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SUMMARY 


A  numerical  and  experimental  study  of  time  domain  methods  for  modeling  and 
parameter  identification  of  structural  systems  is  presented.  Models  are  developed  which 
can  be  used  to  predict  the  transient  response  of  multiple  degree  of  freedom  systems 
subjected  to  arbitrary  input.  The  linear,  discrete  time  transfer  function  is  expressed  in  a 
form  called  the  Autoregressive  Moving  Average  (ARMA)  model.  The  ARMA  model  is  a 
minimum  parameter  model  that  may  be  parameterized  with  a  minimum  number  of  measured 
quantities.  The  ARMA  model  is  contrasted  to  traditional  models  such  as  differential 
equation  models  and  modal  methods.  The  ARMA  model  identification  algorithms  are  also 
compared  to  time  domain  modal  methods.  Though  it  has  proven  useful  for  a  number  of 
low  order  systems,  the  identification  of  the  ARMA  model  is  often  hampered  by  the 
sensitivity  of  parameter  estimates  to  measurement  noise  bias.  Modal  methods  and  signal 
processing  algorithms  reduce  the  sensitivity  to  noise  by  using  methods  such  as  model 
overspecification  or  complex  iterative  procedures.  Numerous  algorithms  and  their 
effectiveness  in  reducing  noise  bias  are  presented.  Finally,  a  two-stage  identification 
technique  is  presented  which  uses  overspecified  models  to  estimate  the  poles  of  the  ARMA 
model  from  system  free  vibration  response  data  using  a  form  of  the  Principal  Eigenvectors 
Method.  A  Least  Squares  algorithm  is  then  used  to  identify  the  transfer  function  zeroes 
from  forced  response  records.  Experimental  response  data  for  both  single-degree-of- 
freedom  and  multiple-degree-of-freedom  systems  are  used  to  evaluate  the  different 
methods  of  identification.  An  extension  of  the  methods  to  nonlinear  discrete  time  series 
models  is  al.so  presented.  These  models  are  a  logical  extension  of  the  linear  models.  The 
nonlinear  model  forms  are  discussed  and  both  simulated  and  actual  experimental  data  from 
an  oleo-pneumatic  landing  gear  strut  are  used  to  evaluate  the  model  forms. 
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INTRODUCTION 


Background 

Dynamic  system  modeling  of  real  structural  systems  is  typically  done  for  one  of 
three  reasons.  The  system  model  may  be  needed  to  determine  whether  the  structure  will 
have  acceptable  response  to  the  variety  of  input  forces  that  it  will  encounter  during  its 
lifetime.  If  the  response  is  unacceptable,  either  the  structure  is  modified,  perhaps  using 
some  information  extracted  from  the  model,  or  the  input  forces  have  to  be  altered.  The 
second  reason  is  that  a  system  model  may  be  needed  to  implement  a  control  algorithm.  The 
required  accuracy  of  the  model  depends  on  the  sophistication  of  the  algorithm.  The  third 
reason  may  be  simply  a  verification  of  results  developed  from  some  other  analysis. 

Analytical  or  experimental  procedures  are  used  to  develop  the  system  model.  The 
traditional  analytical  procedures  involve  the  formulation  of  differential  equations  termed  the 
equations  of  motion.  Either  numerical  or  analytical  integration  is  used  to  determine  the 
solutions  of  the  differential  equations  of  motion  and  thus  the  system  response.  One  of  the 
more  popular  approaches  to  develop  the  equations  of  motion  for  a  structure  is  the  Finite 
Element  Method  (FEM).  The  FEM  can  be  used  to  develop  a  set  of  second  order  differential 
equations  through  discretization  of  the  structure.  The  set  of  second  order  differential 
equations  represents  one  form  of  a  system  model.  Structural  engineers  usually  prefer  the 
system  model  in  an  alternative  form  termed  the  modal  form.  The  modal  form  of  the  system 
model  is  comprised  of  the  mode  shapes,  modal  frequencies,  and  modal  damping  factors. 
The  modal  frequencies  and  modal  damping  factors  can  be  determined  from  the 
characteristic  equation  of  the  structure.  The  mode  shape  defines  the  amplitude  of  response 
at  given  modal  frequencies  at  the  discretized  points.  The  response  at  the  discrete  point  is  a 
summation  of  the  scaled  response  at  each  modal  frequency.  The  modal  form  of  the  model 
can  be  synthesized  from  the  equations  of  motion  and  the  equations  of  motion  can  be 
synthesized  from  the  modal  form  of  the  model. 

The  traditional  experimental  procedure  for  the  determination  of  a  system  model 
involves  the  modal  form.  The  modal  information  is  derived  from  experimentally  acquired 
data.  The  structure  must  be  excited  at  a  point  or  a  number  of  points  and  the  response 
measured  at  a  discrete  set  of  points.  Frequency  domain  methods  of  identification  transform 
the  time  domain  response  measurements  into  the  frequency  domain  in  the  form  of  Discrete 
Fourier  Transforms  (DFT)  using  a  Fast  Fourier  Transform  (FFT)  algorithm.  Frequency 
Response  Functions  (FRF)  can  be  derived  by  dividing  the  DFT's  of  the  response 
measurements  by  the  DFT's  of  the  input  functions.  Multiple  data  sets  are  acquired  and  the 
FRF's  for  each  input-output  pair  are  averaged  to  reduced  leakage  and  measurement  errors. 
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The  modal  information  is  then  parameterized  by  curve  fitting  techniques  performed  on  the 
ensemble  of  FRF's. 

In  the  past  decade,  modal  analysis  techniques  using  time  domain  data  have  also 
become  piopular.  These  techniques  use  time  domain  data  to  parameterized  the  modal  form 
of  the  model.  Either  a  set  of  time  domain  free  response  data  can  be  directly  used  or  the 
ensemble  of  FRF's  can  be  transformed  back  into  the  time  domain  providing  impulse 
response  functions.  The  major  time  domain  methods  include  the  Complex  Exponential 
Method  (CEM),  Ibrahim's  Time  Domain  technique  (ITD),  and  the  Eigensystem  Realization 
Algorithm  (ERA).  The  time  domain  methods  involve  the  formulation  of  a  coefficient 
matrix  from  the  time  domain  data.  The  eigenvalues  and  eigenvectors  of  this  matrix 
completely  characterize  the  modal  form  of  the  model.  The  methods  differ  by  the  way  in 
which  the  coefficient  matrix  is  formed. 

Time  domain  methods  offer  a  higher  resolution  capability  for  closely  spaced 
frequencies  than  the  frequency  domain  methods.  The  major  advantage  of  the  frequency 
domain  methods  is  that  multiple  sets  of  data  can  be  used  and  the  effect  of  noise  can  be 
greatly  reduced.  The  time  domain  methods  have  to  use  overspecified  model  orders  or 
iterative  techniques  to  counter  the  effects  of  noise  bias.  The  frequency  domain  methods  are 
user  intensive,  requiring  the  user  to  distinguish  between  important  and  unimportant  modes; 
the  time  domain  methods  are  more  automated.  The  time  domain  methods  may  therefore  be 
better  suited  for  rapid  "on-line"  modeling. 

The  frequency  domain  models  can  be  represented  by  a  matrix  of  transfer  functions 
in  the  Laplace  domain.  The  time  domain  models  can  be  represented  by  a  matrix  of  discrete 
time  transfer  functions.  Response  is  determined,  using  the  frequency  domain  models,  by 
convolution  integrals  involving  the  elements  of  the  matrix  and  the  input  function. 
Response  is  determined,  using  the  time  domain  models,  by  recursive  addition  of  state 
variables.  The  weights  of  the  state  variables  in  the  addition  process  are  defined  by  the 
parameters  in  the  discrete  time  transfer  function. 

The  Problem 

During  service,  certain  types  of  structures  can  assume  many  different 
"configurations"  due  to  changes  in  mass  loading,  temperature  changing  the  structural 
properties,  or  failure  of  localized  components.  If  one  is  concerned  with  the  prediction  of 
the  dynamic  response  of  critical  components  of  a  complete  structure  when  subjected  to  a 
specific  loading  condition,  then  a  single  structural  model  defined  for  the  "ideal"  structure 
may  not  prove  suitable.  The  many  possible  structural  "configurations"  and  inputs  in  the 
lifetime  of  some  structures  place  limits  on  the  usefulness  of  certain  detailed  analytic  or 
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numerical  modeling  procedures  such  as  Finite  Element  Analysis.  It  would  be  desirable  to 
have  techniques  that  would  provide  "on-line"  or  adaptive  models  that  would  be  suitable  for 
the  rapid  and  accurate  prediction  of  transient  dynamic  response  for  a  specific  set  of  loading 
conditions. 

Consider  the  problem  of  an  airplane  about  to  traverse  a  rough  or  damaged  runway. 
Prediction  of  peak  acceleration  levels  at  several  critical  locations  on  the  stmcture  could  help 
to  determine  whether  the  aircraft  would  be  damaged  while  using  the  runway  or  could  help 
to  determine  the  speed  and  path  the  aircraft  should  take  to  minimize  the  risk  of  dsimage.  A 
model  relating  the  vibration  (peak  displacement  or  acceleration)  at  the  critical  locations 
within  the  structure  to  the  specific  runway  profile  would  be  needed.  However,  each  time 
the  aircraft  is  used,  its  structural  configuration  may  change;  for  example,  the  payload  or 
fuel  loading  may  be  modified.  When  the  structural  configuration  changes,  the  predictive 
model  must  change.  Or  consider  future  large  space  structures  that  will  be  fabricated  in 
orbit.  Dynamic  models  will  be  necessary  to  control  these  structures.  The  varying  structural 
characteristics  during  assembly  would  outdate  the  baseline  analytical  or  modal  models  used 
to  predict  dynamic  response. 

In  both  examples,  input  to  the  structure  is  provided  at  a  discrete  number  of  points, 
and  knowledge  of  the  system  response  may  be  necessary  at  only  a  finite  number  of  points 
on  the  structure.  It  would  be  desirable  to  identify  or  update  the  structural  characteristics 
every  time  the  structure  is  used.  Ideally,  digital  time  histories  of  the  input  and  system 
output  would  be  acquired,  and  the  system  model  or  transfer  function  would  be  derived  or 
updated  automatically.  For  linear  systems,  discrete  time  transfer  functions  are  well-suited 
for  automatic  digital  identification.  One  form  of  a  discrete  time  transfer  function  is  the 
Autoregressive  Moving  Average  (ARMA)  model.  The  ARMA  model  can  be  used  to  predict 
response  for  Multiple  Input  Multiple  Output  (MIMO)  systems  as  well  as  Single  Input 
Single  Output  (SISO)  systems. 

Overview 

This  report  is  concerned  with  the  suitability  and  identification  of  ARMA  models  for 
structural  response  prediction,  and  will  be  limited  to  SISO  applications.  The  extension  to 
MIMO  systems  seems  to  be  direct.  (More  research  in  this  area  may  be  necessary.)  This 
report  presents  a  discussion  of  the  suitability  of  using  ARMA  models  for  structural 
response  prediction  and  an  overview  of  the  factors  relating  to  the  parameter  identification 
process.  The  ARMA  model  was  chosen  because  of  its  suitability  to  the  digital  environment 
and  because  time  domain  models  can  be  identified  in  a  more  automated  fashion  than 
frequency  domain  models.  These  two  factors  are  very  important  when  considering  the 
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ultimate  goal  of  response  prediction.  For  example,  if  an  aircraft  nxi  response  is  desired, 
the  data  required  to  parameterize  a  model  for  the  aircraft  would  be  acquired  as  part  of  the 
digital  instrumentation  system  planned  for  future  aircraft  systems,  and  the  "model"  for 
structural  response  prediction  would  have  to  be  automatically  synthesized  so  that 
"warnings"  or  operational  instructions  would  be  provided  to  the  pilot. 

In  addition  to  the  ARMA  model.  Autoregressive  (AR)  models  and  Moving  Average 
(MA)  models  will  be  discussed.  The  major  identification  algorithms  for  the 
parameterization  of  the  AR,  MA,  and  ARMA  models  will  be  presented.  A  classical 
identification  scheme,  the  single-stage-least-squares  (SLS)  algorithm,  and  a  proposed 
modification,  a  two-stage-least-squares  (2LS)  algorithm,  are  the  main  algorithms 
presented.  Results  using  the  SLS  and  2LS  algorithms  to  parameterize  single-degree-of- 
freedom  (SDOF)  and  multiple-degree-of-freedom  (MDOF)  system  models  are  presented. 
These  results  were  achieved  using  simulated  and  actual  experimental  data  and  the  models 
were  then  used  to  predict  transient  response  for  these  systems.  The  final  section  will 
describe  the  extension  of  the  discrete  time  series  models  to  nonlinear  models. 
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DYNAMIC  SYSTEM  MODELS 


There  are  many  different  types  of  analytic  or  numerical  models  that  can  be  used  to 
predict  vibration.  The  models  can  be  developed  from  either  basic  physical  principles  or 
from  the  observed  response  of  the  system.  This  report  will  focus  on  the  development  of 
numerical  models  for  complex  structural  systems  by  recording  and  analyzing  the  response 
of  the  system.  The  experimental  model  can  be  in  the  form  of  a  parameterized  set  of 
differential  equations,  a  modal  model,  or  a  parameterized  set  of  difference  equations.  This 
section  briefly  presents  the  different  forms  of  the  system  model  and  compares  the 
requirements  necessary  to  determine  system  parameters.  A  discussion  at  the  end  of  the 
section  will  give  some  of  the  reasons  why  the  ARMA  model  has  been  chosen  as  the  focus 
of  this  report. 

Differential  Equations 

The  undergraduate  engineer  is  taught  in  mechanics  and  vibration  classes  that 
complex  structures  can  be  modeled  as  discretized,  linearized,  multiple-degree-of-ffeedom 
systems.  Newtonian  mechanics  can  be  used  to  write  a  system  of  differential  equations  to 
describe  the  dynamic  behavior  of  the  structure.  The  differential  equations  are  typically 
written  in  the  form  of  second  order,  linear,  constant  coefficient  equations  as 

M  y  +  Cy  +  Ky  =  Bu  (1) 

where  y  is  a  vector  representing  the  position  of  n  generalized  discrete  points  representing 
the  n  degrees  of  freedom  and  u  is  a  vector  of  m  inputs.  The  coefficient  matrices  M,  C,  K 
and  B,  which  are  referred  to  as  system  parameters,  are  used  to  describe  the  system.  Some 
parameter  identification  techniques  attempt  to  determine  M,  C,  K,  and  B  from 
experimental  data.  The  estimation  procedures,  referred  to  as  equation  of  motion  methods, 
require  that  all  inputs  and  the  position,  velocity,  and  acceleration  at  the  generalized  degrees 
of  freedom  be  measured.  Some  methods  use  integration  or  differentiation  to  generate  two 
of  three  measurement  quantities  (position,  velocity,  and  acceleration),  from  one 
measurement  (position,  velocity,  or  acceleration)  but  the  applicability  of  this  approach  to 
anything  but  the  simplest  of  systems  has  not  been  demonstrated.  Usually,  the  matrices  in 
Eqn  1  can  only  be  estimated  to  within  a  transformation  unless  a  mass  or  stiffness 
perturbation  method  is  used.  If  the  parameters  of  the  system  can  be  determined  and  one 
wants  to  determine  the  response  of  the  system  to  another  input,  the  usual  methods  for  the 
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solution  of  linear  differential  equations  are  sufficient.  Given  the  model,  the  response  of  the 
generalized  points  io  an  arbitrary  input  can  be  predicted  through  numerical  integration. 

The  set  of  differential  equations  can  be  reduced  to  a  single  differential  equation  of 
order  2n  with  a  single  dependent  response  variable.  The  dependent  variable  can  be  any 
linear  combination  of  the  n  generalized  degrees  of  freedom.  The  single  differential  equation 
can  also  be  parameterized  from  experimental  data.  The  estimation  procedure  requires  that 
all  inputs  and  the  dependent  response  with  all  of  its  derivatives  up  to  order  2n  be  measured 
or  estimated.  The  differential  equation  can  be  used  to  predicted  the  response  of  the  single 
dependent  variable  to  an  arbitrary  input.  The  reduction  of  the  model  given  by  Eqn  1 
simplifies  the  estimation  process,  since  the  number  of  parameters  of  the  model  has  been 
reduced.  The  reduced  model  gives  response  predictions  of  the  one  dependent  variable  and 
not  at  all  the  generalized  coordinates. 

Modal  Models 

Modal  methods  use  the  measured  response  to  estimate  vibration  parameters  such  as 
the  mode  shape,  natural  frequencies,  damping  factors,  and  modal  participation  factors.  The 
input  to  the  system  can  take  the  form  of  an  impulse,  a  "random"  time  history,  or  a  periodic 
function  of  specific  form,  such  as  a  sine  wave.  The  modal  method  assumes  that  the  actual 
system  response  is  the  superposition  of  a  number  of  characteristic  "modes"  with  specific 
frequencies  and  damping  characteristics.  The  method  is  based  on  the  measurement  of  the 
complex  response  of  the  system  and  then  the  processing  of  the  response  data  to  determine 
the  contributing  modal  parameters. 

The  response  of  the  system  due  to  an  impulse  is  the  free  response  after  the  impulse 
has  occurred.  The  solution  of  Eqn  1  due  to  an  impulse  is  of  the  form  y(t)=p  exp(>.t). 
Substitution  into  Eqn  1  yields  the  eigenvalue  problem 

[X2  M  +  X  C  +K]  p=  0  (2) 

There  exist  2n  discrete  time  eigenvalues  (poles),  X,  representing  the  natural  frequencies  and 
damping  factors.  Here,  p  represents  the  complex  modes  of  vibration  and  is  a  scaled  version 
of  the  mode  shapes.  The  system  response  to  the  unit  impulse  is  found  by  the  summation 

2n 

X(0=X  Pj 

j=l 
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Modal  methods  use  time  domain  data  or  frequency  domain  data  to  identify  the  poles  and  the 
complex  modes  of  vibration.  The  response  at  the  measurement  points  due  to  an  impulse  at 
the  other  input  locations  is  needed  to  identify  the  complex  modes  of  vibration  for  each 
input;  the  scaling  of  the  individual  modes  define  the  modal  participation  factors. 

The  identification  of  the  modal  model  requires  that  the  impulse  response  function, 
in  the  time  domain,  or  the  frequency  response  function,  in  the  frequency  domain,  be 
"measured"  at  every  output  location  for  every  input  location.  Whether  the  identification 
procedure  is  time  or  frequency  based,  the  impulse  response  function  is  usually  developed 
in  the  frequency  domain  using  ensemble  averages  of  the  input  and  output  frequency 
spectra.  A  measurement  of  the  response  at  ever>'  location  and  a  measurement  of  each  input 
are  required.  The  measurements  can  be  made  all  at  once  or  separately.  But  many 
measurements  at  many  locations  are  usually  needed.  This  method  is  very  well-established 
and  is  very  useful  if  a  complete  model  of  the  structure  is  desired. 

There  are  some  advantages  and  disadvantages  to  the  use  of  modal  methods.  The 
estimated  frequencies,  damping  factors,  and  mode  shapes  from  each  impulse  response  will 
be  slightly  different.  Global  estimates  have  to  be  made  through  the  use  of  some  averaging 
procedure.  Also,  if  there  are  n  degrees  of  freedom,  response  measurements  at  n  locations 
are  required.  Given  frequency,  damping  factor,  mode  shape  estimates,  and  modal 
participation  factors,  the  generalized  stiffness  (M'^K),  generalized  damping  matrices  (M" 
^C),  and  generalized  input  matrix  can  be  estimated.  Identification  of  the  full 

modal  model  is  equivalent  to  the  identification  of  the  differential  equations  when  the 
generalized  coordinates  are  the  modes  of  vibration.  The  modal  model  can  be  synthesized 
from  the  differential  equation  model  and  vice  versa. 

Difference  Models 

The  differential  equation  model  can  be  transformed  into  a  difference  equation  model 
by  a  number  of  differencing  techniques.  A  form  of  the  difference  equation  model  is 

-»■  A,  y^.i  +  A2  yk.2=  Bp  u,  Bp  2 

where  y  is  a  vector  of  n  responses  and  u  is  a  vector  of  m  inputs  at  discrete  times.  Equation 
4  is  an  ARMA  model  for  an  Multiple  Input  and  Multiple  Output  (MIMO)  system.  The 
Autoregressive  parameters  are  the  A  matrices  and  Moving  Average  parameters  are  the  B 
matrices.  Mode  shapes,  frequencies,  and  damping  factors  can  all  be  derived  from  a  MIMO 
ARMA  model  as  described  by  reference  1.  Identification  of  the  ARMA  model  given  by 
Eqn  4  would  require  the  input  at  each  location  to  be  measured  and  the  response  at  n 
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locations  to  also  be  measured,  just  as  with  the  modal  models.  But  unlike  the  modal 
models,  the  ARMA  model  can  be  identified  without  forming  the  impulse  response 
functions.  Reference  1  also  shows  that  identification  of  the  Autoregressive  matrices  from 
impulse  response  functions  is  equivalent  to  the  time  domain  modal  method  known  as  the 
Ibrahim  Time  Domain  Technique  (reference  2).  Unlike  the  differential  equation  models,  no 
integration  is  necessary  to  predict  response.  But  like  the  differential  equation  models,  the 
ARMA  model  can  be  reduced  if  only  a  single  response  is  measured.  The  resulting 
difference  model  would  be 

y,  +  a,y^_j+...+a2^y,_^^=b/  u,+  hj  (6) 

which  is  a  Multiple  Input  Single  Output  (MISO)  model;  here  the  Moving  Average 
parameters  are  row  vectors  and  the  Autoregressive  parameters  are  constants.  For  a  single 
input,  the  Moving  Average  parameters  become  scalars.  The  identification  of  the  SISO 
ARMA  model  requires  that  only  a  single  input  and  a  single  output  be  measured. 

Overview 

The  goal  of  the  research  effort  summarized  in  this  report  was  to  study  methods  that 
may  be  able  to  automatically  parameterize  a  model  for  a  complex  structural  system.  The 
model  would  then  be  used  to  predict  response  for  a  specified  system  input.  The 
Differential  Equation  methods  are  undesirable  because  the  full  state  of  the  response 
(position, velocity,  and  acceleration)  is  required  from  measurements  or  from  auxiliary 
estimations  (integrations  or  differentiations).  The  modal  methods  are  undesirable  because 
multiple  measurement  locations  are  required.  The  ARMA  model  appears  to  present  the 
best  alternative.  The  size  of  the  ARMA  models  is  adjustable,  depending  upon  how  many 
response  locations  are  required.  In  the  aircraft  taxi  response  example,  if  the  predicted 
acceleration  response  at  the  wingtip  and  the  pilot  location  are  the  only  critical  locations,  the 
ARMA  model  can  be  developed  for  the  response  at  just  those  locations.  Only  the 
acceleration  at  wingtip  and  pilot  locations  and  the  inputs  have  to  be  measured  to  identify  the 
model.  The  ARMA  model  can  also  be  identified  without  forming  the  impulse  response  in 
the  frequency  domain.  The  identification  of  the  ARMA  model  can  be  accomplished  entirely 
in  the  time  domain.  The  next  section  will  describe  the  ARMA  model  in  depth,  and  the 
remaining  sections  deal  with  the  identification  of  the  discrete  time  series  models. 
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DISCRETE  TIME  TRANSFER  FUNCTIONS 


The  ARMA  model  is  used  in  this  report  to  predict  structural  vibration. 
Identification  procedures,  examined  later  in  the  report,  will  parameterize  the  ARMA  model 
from  input-output  discrete  time  histories.  Before  proceeding  to  the  identification 
procedures,  background  on  the  ARMA  model  is  provided.  This  section  describes,  through 
the  use  of  the  z-transform,  the  discrete  time  transfer  function  and  its  relationship  to  the 
continuous  time  transfer  function.  Approximations  are  introduced  that  will  allow  the 
discrete  time  transfer  function  to  be  synthesized  from  the  continuous  time  transfer  function. 
Applications  of  the  ARMA  model  are  briefly  discussed,  and  the  advantages  and 
disadvantages  of  ARMA  models  are  presented.  Examples  using  discrete  time  transfer 
functions  to  predict  vibration  are  included.  This  section  introduces  the  ARMA  model, 
describes  relationships  to  continuous  transfer  functions,  and  shows,  by  example,  that  the 
ARMA  model  can  be  used  to  accurately  predict  structural  vibration. 

The  Z-transform 

For  linear  systems  the  relationship  between  the  input  to  the  system  and  the  system 
response  or  output  can  be  expressed  as  transfer  functions.  The  transfer  functions  can  be 
determined  by  transforming  the  input  and  output  time  functions  using  Laplace 
transformations.  The  transfer  function  in  the  Laplace  domain  is  found  by  dividing  the 
Laplace  transform  of  the  output  by  the  Laplace  transform  of  the  input, 

,  Y(s) 

“  U(s) 

The  frequency  domain  transfer  function  is  simply  the  Laplace  domain  transfer  function 
evaluated  along  the  complex  axis  (s=ja)).  The  goal  of  frequency  domain  identification 
methods  is  to  accurately  determine  a  frequency  domain  transfer  function  for  every  input- 
output  pair.  There  exists  a  discrete  time  equivalent  of  Eqn  7  that  can  be  used  in  the  time 
domain.  The  discrete  time  equivalent  involves  the  use  of  the  z-transform  instead  of  the 
Laplace  transform.  The  relationship  between  the  s  or  Laplace  domain  and  the  z  or  discrete 
time  domain  is  written  as 

z=e^*^  (8) 

where  s  is  the  Laplace  variable,  z  is  the  z-transform  variable,  and  h  the  disc’-ete  time 
interval.  Note  that  the  stable  region  in  the  Laplace  domain  ( that  is  the  region  in  which  the 
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real  part  of  s  in  negative)  is  mapped  inside  the  unit  circle  in  the  z-domain.  Analytically,  the 
discrete  time  transfer  function  is  obtained  by  a  straightforward  procedure  if  the  input  and 
output  time  functions  are  known.  Realistically,  the  time  functions  are  measured  at  discrete 
time  instances,  and  the  exact  input  and  output  time  functions  are  unknown.  In  practice,  the 
discrete  time  transfer  function  has  to  be  determined  using  a  modified  form  of  the  z- 
transform. 

Suppose  that  a  function,  f(t),  is  known  precisely  at  n  discrete  time  instances,  then 
the  sampled  function,  f* (t),  can  be  reconstructed  as  a  summation  of  the  discrete  time  values 
times  the  appropriate  delta  functions, 

f^(t)  =  ^  f(kh)  5(t  -  kh)  (9) 

k=0 


Here,  f* (t)  is  exactly  equivalent  to  f(t)  at  the  n  discrete  time  instances  and  is  zero  at  every 
other  time.  The  Laplace  transform  of  Eqn  9  is 

F(s)  =  ^  f(kh)  e  '“‘’  (10) 


The  z-transform  of  Eqn  9  is  calculated  by  substituting  Eqn  8  into  Eqn  10  to  obtain 

F(z)  =  ^  f(kh)  z (11) 

k=0 


The  discrete  time  transfer  function  is  simply  the  ratio  between  the  z-transform  of  the  output 
and  the  z-transform  of  the  input; 


n-1 


G(z)  = 


^  y(kh)  z'*^ 

Y(z)  fro 
U(z)"  n-1 


(12) 


If  there  is  no  noise  or  error  in  the  measurements,  Eqn  12  may  reduce  to  a  polynomial 
expression 
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,  U  -I  u  '2  u  -p 

+  bjZ  +  h2^  +  ...+  bpZ  ^ 

G(z)  - - ^ - :5 - — 

1  +  ajz  +  +  ...+  apZ  ^ 


(13) 


If  there  is  noise,  an  identification  algorithm  must  be  used  to  estimate  the  parameters  in  Eqn 
13.  The  discrete  time  transfer  function  shown  in  Eqn  13  is  also  called  the  Infinite  Impulse 
Response  (HR)  filter  and  can  be  written  in  the  form  of  an  ARMA  model.  The  basic 
ARMA  model  is  the  SISO  model  which  has  the  form 


y(k)=  -aiy(k-l)  -a2y(k-2)  -...-apy(k-p)  +b()u(k)  +biu(k-l)  +...+bpu(k-p)  (14) 

where  u(k)  are  the  di.screte  values  (at  time  kh)  of  the  various  input  functions,  and  y(k)  are 
the  di.screte  time  values  of  the  response.  The  variable  z  in  Eqn  13  is  also  known  as  the 
forward-shift  operator.  The  forward  shift-operator  shifts  the  discrete  time  index  by  1 , 

z  [y(k)]  =  y(k+l)  (15) 

The  discrete  lime  transfer  funcMon  is  only  an  approximation  to  the  continuous  time  transfer 
function.  Variations  in  the  input  between  sample  instances  do  not  affect  the  output 
according  to  the  discrete  time  transfer  function,  although  such  variations  will  affect  the 
output  of  the  actual  system.  Notice  that  two  different  time  functions  can  have  exactly  the 
same  z-transform. 

An  alternative  formulation  of  the  discrete  time  transfer  function  is  given  in 
Appendix  A.  The  derivation  begins  with  the  state  space  representation  of  a  system  and 
ends  with  the  discrete  time  transfer  function.  The  only  approximation  used  in  the 
derivation  is  the  assumption  that  the  input  is  constant  over  the  sample  interval. 

Transformation  to  the  Discrete  time  domain 

Ideally,  to  transform  a  transfer  function  in  the  Laplace  domain  to  the  time  domain, 
the  substitution 

s  =  '^lnz  (16) 

would  conven  the  Laplace  transforms  to  z-transforms.  The  main  difficulty  with  the 
conversion  is  that  the  transfer  function  will  not  reduce  to  a  simple  form.  The  resulting  form 
will  contain  terms  involving  the  logarithm  of  z,  which  do  not  have  the  convenient  definition 
of  the  forward-shift  operator,  Eqn  7.  One  would  like  to  have  a  substitution  approximately 
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equivalent  to  Eqn  16  that  would  result  in  discrete  time  transfer  functions  that  were  ratios  of 
polynomials  in  z.  Such  substitutions  are  commonly  used  in  digital  control  engineering. 
The  simplest  are  the  forward  rectangular  rule. 


and  the  backward  rectangular  rule. 


The  most  popular  approximation  is  Tustin's  rule. 


^"hz+1 


The  forward  rectangular  rule  is  the  Euler  integration  scheme.  TTie  backward  rectangular  rule 
corresponds  to  backward  differencing  and  Tustin's  rule  is  trapezoidal  or  Modified-Euler 
integration.  The  input  is  assumed  to  be  constant  in  the  forward  and  backward  rectangular 
rules  and  varies  linearly  in  Tustin’s  rule.  The  rules  are  approximations  of  Eqn  8,  the 
expansion  of  Eqn  8  is 


sh 

z  =  e  = 


I 


•  =  1  +  sh  + ' 


(sh)^  (sh)^ 


The  forward  rectangular  rule  is  the  two-term  approximation  of  Eqn  8.  The  backward 
rectangular  rule  is  expressed  as 


which  is  a  two-term  approximation  plus  an  over  approximation  of  the  remaining  terms. 
Tustin's  rule  is  expanded  as 


'4 

z=.^=l+sh  + 

'■T 


which  is  a  three-term  approximation  plus  an  over  approximation  of  the  remaining  terms. 

The  reason  that  Tustin’s  rule  is  the  most  popular  approximation  is  that  it  maps 
stable  poles  in  the  Laplace  domain  to  stable  pole  locations  in  the  z-domain  and  it  maps 


unstable  poles  in  the  Laplace  domain  to  unstable  pole  locations  in  the  z-domain.  The 
forward  rectangular  rule  maps  unstable  poles  in  the  Laplace  domain  to  unstable  pole 
locations  in  the  z-domain  but  maps  some  stable  poles  in  the  Laplace  domain  to  unstable 
pole  locations  in  the  z-domain.  Tlie  backward  rectangular  rule  maps  stable  poles  in  the 
Laplace  domain  to  stable  pole  locations  in  the  z-domain  but  maps  some  unstable  poles  in 
the  Laplace  domain  to  stable  pole  locations  in  the  z-domain.  References  3  and  4  describe 
these  approximations  in  some  detail. 

Applications  of  ARMA  models 

The  accuracy  of  the  approximations  discussed  in  the  previous  section  improves  as 
the  sample  rate  increases.  All  three  provide  nearly  the  same  discrete  time  transfer 
functions.  The  approximations  have  received  extensive  use  in  digital  control  engineering. 
A  crude  procedure  in  control  engineering  is  to  analytically  design  a  compensator  in  the 
Laplace  domain,  transform  it  at  a  high  sampling  rate  to  the  z-domain,  and  then  use  a 
controller  to  implement  the  discrete  time  transfer  function  of  the  compensator  (reference  4). 

The  primary  application  of  the  methods  discussed  in  this  report  is  not  system 
control  but  is  focussed  on  the  discrete  time  transfer  functions  or  ARMA  models  when  used 
for  vibratory  prediction  of  complex  multiple-degree-of-freedom  systems.  The 
approximations  are  useful  in  discovering  what  considerations  are  necessary  to  make  the 
discrete  time  transfer  function  an  attractive  and  accurate  means  of  predicting  vibration.  The 
approximations  suggest  that  a  high  sample  rate  may  be  needed  to  increase  the  accuracy, 
although,  a  high  sample  rate  increases  the  computational  load  of  the  prediction. 

Other  control-oriented  applications  of  the  ARMA  model  include  digital  filtering 
(references  5  and  6)  and  identification  of  chemical  processes  (reference  7).  Structural 
applications  of  the  ARMA  model  include  substnicture  modelling  (reference  8),  probabilistic 
simulations  (reference  9),  spectral  estimation  (reference  10),  and  vibratory  parameter 
extraction  (reference  1 1). 

Advantages  and  Disadvantages  of  the  ARMA  model 

There  are  advantages  to  the  u.se  of  the  ARMA  models  in  identification  and 
prediction.  The  ARMA  model  is  a  minimum  parameter  model;  it  therefore  requires  the  least 
number  of  parameters  to  characterize  the  system.  Consider  an  N  degree  of  freedom,  linear 
dynamic  system  which  can  be  modeled  as  N  second  order,  ordinary  differential  equations. 
The  2N  states  of  this  model  can  be  expressed  as  the  generalized  coordinates  and  their  first 
time  derivatives.  A  fully  coupled  system  with  viscous  damping  would  require  3(N2-(-NV2 
parameters  (assuming  symmetric  mass,  damping  and  stiffness  matrices)  to  completely 
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define  the  system.  Solution  of  this  system  of  equations  would  provide  the  complete  spatial 
distribution  for  the  structure  at  any  instant  in  time.  An  alternate  form  of  the  system  model 
could  be  a  single  ordinary  differential  equation  of  order  2N.  The  states  of  this  model  would 
be  a  single  generalized  coordinate  and  all  of  its  derivatives  up  to  order  2N.  For  free 
vibration  response,  this  equation  could  be  characterized  using  only  2N-1  parameters.  The 
solution  would  yield  a  single  state  and  all  of  its  derivatives  up  to  order  2N.  The  second 
form  is  a  minimum  parameter  model.  An  equivalent  ARMA  model  would  also  be  a 
minimum  parameter  model. 

The  SISO  ARMA  model's  states  are  discrete  time  values  of  a  single  input  and  single 
output  at  a  given  instance  in  time  and  at  a  finite  number  of  previous  instances.  The  model 
state  for  a  pth  order  model  contains  the  system  output  at  p  previous  time  instances  and  the 
system  input  at  p+1  time  instances.  The  model  state  at  the  next  time  instance  is  found  by 
incrementing  all  the  states  at  the  current  time  by  1.  The  identification  and  prediction 
procedures  for  the  ARMA  model  require  the  minimum  number  of  measured  quantities. 
Only  the  system  input  and  response  need  be  measured.  For  instance,  if  the  response  is 
acceleration  at  some  location,  velocity  at  that  location  need  not  be  known  to  identify  the 
model.  Once  the  parameters  of  the  model  are  known,  prediction  can  be  quickly 
accomplished,  given  the  input  time  history  and  the  initial  state  of  the  model.  The 
disadvantages  of  ARMA  models  are  their  limited  accuracy,  it  is  an  approximation,  and  the 
identification  problems  associated  with  time  domain  procedures. 

Examples 

In  the  following  examples,  a  single-degree-of-freedom  (SDOF)  base  motion 
problem  and  an  arbitrary  multiple-degree-of-freedom  (MDOF)  system  are  used  to  examine 
the  effect  of  sample  rate  upon  the  accuracy  of  the  discrete  time  transfer  function.  Consider 
the  SDOF  spring-mass-damper  base  motion  problem  shown  in  Figure  1,  the  input  is 
provided  by  the  base  mass  and  the  output  is  the  response  of  the  free  mass.  The  transfer 
function  of  this  simple  linear  system  is 

2Ctos  to^ 

G(s)  =  - 5-  (23) 

S  +2^0)5  +  a 

A  discrete  time  transfer  function  can  be  found  through  any  of  the  approximations  di.scussed 
above.  Tustin's  rule  (Eqn  19)  will  be  used  for  the  example  presented  here.  The  parameters 
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and  accuracy  of  the  discrete  time  transfer  function  depend  upon  the  sample  rate.  It  is 
convenient  to  compare  the  different  transfer  functions  through  the  use  of  the  magnitude 
frequency  response  plot.  The  frequency  response  is  found  by  the  substitution  of  s=j(27if) 
into  the  transfer  function.  The  magnitude  frequency  response  plot  is  found  by  plotting  the 
magnitude  of  resulting  complex  expression  versus  frequency,  f.  The  magnitude  frequency 
response  plot  shows  graphically  the  output  gain  of  the  system  versus  the  input  frequency. 
The  phase  of  the  frequency  response  also  contains  valuable  information,  but  the  meaning 
of  the  magnitude  plot  is  easier  to  interpret  and  will  be  used  for  illustrative  purposes.  Note 
that  the  magnitude  frequency  response  plot  of  the  discrete  time  transfer  function  requires 
the  substitution  of  Eqn  2  prior  to  the  substitution  of  s=j(27tf).  The  transformation  from  the 
discrete  time  domain  to  the  s  domain  is  not  an  approximation.  Also  the  frequency  response 
function  for  the  discrete  time  equivalent  can  be  evaluated  at  any  frequency  value  by  the 
computer. 

The  natural  frequency  of  10  Hz  and  damping  factor  of  1  percent  critical  were  used 
in  Eqn  23  for  the  example.  The  magnitude  frequency  response  plot  of  the  transfer  function 
is  shown  in  Figure  2  along  with  the  magnitude  frequency  response  plots  of  the  discrete 
time  equivalents  obtained  through  the  use  of  Tustin's  rule  at  two  different  sample  rates. 
The  discrete  time  equivalents  using  the  higher  sample  rate  (200  Hz)  is  approximately 
equivalent  to  the  transfer  function  from  0  to  20  Hz.  The  discrete  time  equivalent  using  the 
lower  sample  rate  (50  Hz)  is  approximately  equivalent  to  the  transfer  function  from  0  to  5 
Hz.  This  would  imply  that  the  sample  rate  of  the  discrete  time  approximation  would  have 
to  be  at  least  10  to  20  times  the  base  bandwidth  of  input  signal  to  get  accurate  vibratory 
predictions. 

An  arbitrary  MDOF  system  was  also  used  to  compare  the  effect  of  sampling  rates 
upon  the  accuracy  of  the  discrete  time  model  for  a  more  complex  system.  The  arbitrary 
system  was  generated  by  randomly  placing  poles  and  zeros  of  the  transfer  function  in  the 
Laplace  domain.  Figure  3  shows  the  magnitude  frequency  response  of  an  arbitrary  system 
(a  4  DOF  system)  along  with  the  magnitude  frequency  response  plots  of  the  discrete  time 
equivalents  obtained  through  the  use  of  Tustin's  rule  at  two  different  sample  rates.  The 
discrete  time  equivalent  at  the  1  kHz  sampling  rate  is  approximately  equivalent  to  the 
arbitrary  transfer  function  over  the  entire  frequency  range  plotted  (0  to  75  Hz).  The 
discrete  time  equivalent  at  the  lower  sampling  rate  (200  Hz)  is  only  accurate  up  to  about  20 
Hz.  Predictions  of  system  response  to  simulated  inputs  using  the  discrete  time  equivalents 
in  the  form  of  ARMA  models  (Eqn  8)  were  compared  to  those  calculated  through  standard 
numerical  integration  (fourth  order  Runge-Kutta).  A  periodic  input  comprised  of  three  unit 
amplitude  sinusoids  at  frequencies  of  2,  30,  and  50  Hz  is  shown  in  Figure  4.  Inspection  of 
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Figure  2:  Tustin’s  Rule  used  to  approximate  a  SDOF  system 

transfer  function. 
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Figure  3:  Tustin's  Rule  used  to  approximate  a  MDOF  system 

transfer  function. 
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Figure  4:  Periodic  input  to  the  MDOF  system. 
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the  magnitude  frequency  response  plots  indicates  that  the  discrete  time  equivalent  sampled 
at  1  kHz  should  have  virtually  the  same  response  as  the  continuous  time  system  at  these 
frequencies  and  the  discrete  time  equivalent  sampled  at  200  Hz  will  not  be  as  accurate  at  the 
higher  frequencies.  Figure  5  shows  the  response  predicted  by  the  ARMA  models  in 
comparison  to  the  numerically  integrated  response.  (The  predictions  are  determined 
through  knowledge  of  the  input  and  model  parameters  only;  the  response  is  calculated 
recursively  and  the  initial  conditions  are  assumed  to  be  zero.)  Figure  5  shows  that  both 
predictions  are  highly  accurate.  One  would  have  expected  the  discrete  time  equivalents  at 
the  higher  sampling  rate  to  be  accurate,  and  the  errors  of  the  discrete  time  equivalent  at  the 
lower  sampling  rate  to  be  small,  since  the  magnitude  of  the  response  is  so  low  at  the  higher 
input  frequencies. 

Figure  6  shows  a  pseudo-random  input  (Gaussian  noise  with  a  unit  variance)  which 
should  provide  a  very  broadband  excitation  for  this  system.  The  bandwidth  of  the  signal 
was  reduced  to  approximately  0-50  Hz  through  the  use  of  a  digital  filter  (an  8th  order 
ARMA  approximation  of  an  8  pole  Butterworth  filter).  The  magnitude  frequency  response 
of  the  filter  is  shown  in  Figure  7.  The  filtered  broadband  input  is  shown  in  Figure  8.  The 
filtered  input  was  then  applied  to  the  s>itcm.  The  predicted  responses,  based  on  the 
discrete  time  equivalents,  are  shown  in  comparison  to  the  numerically  integrated  response 
in  Figure  9.  The  ARMA  approximation  prediction  at  the  high  sampling  rate  follows  the 
numerically  integrated  response  very  closely  (within  the  resolution  of  the  plot).  The 
ARMA  prediction  for  the  lower  sample  rate  is  not  quite  as  accurate,  but  for  many 
applications  it  would  be  considered  to  be  an  adequate  representation  of  the  system 
response. 

Overview 

The  discrete  time  transfer  function  has  been  shown  to  provide  an  accurate  prediction 
of  simulated  response  when  the  sampling  rate  is  at  least  10  to  20  times  larger  than  the 
bandwidth  of  the  input.  The  discrete  time  transfer  functions  used  in  the  examples  were 
developed  through  analytical  approximations  and  the  errors  are  due  to  those 
approximations.  The  discrete  time  transfer  function  assumes  that  the  input  is  constant  over 
the  sampling  interval.  Even  if  a  discrete  time  transfer  function  could  be  found  that  did  not 
use  the  approximations,  the  sampling  rate  would  conceivably  have  to  be  at  least  10  times 
higher  than  the  bandwidth  of  interest  to  realistically  represent  the  input  function.  Although 
Shannon's  Sampling  Theorem  (reference  4)  states  that  only  two  points  per  cycle  are  needed 
to  reconstruct  a  signal,  it  is  well  known  that  the  sampling  rate  needs  to  be  10  to  20  times  the 
bandwidth  upper  limit  to  accurately  perform  numerical  integration  for  most  numerical 
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Figure  5:  Simulated  response  of  the  MDOF  system  to  the  periodic  input. 
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Figure  6:  Pseudo— random  input  to  the  MDOF  system. 
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Figure  8:  Filtered  pseudo-random  input  to  the  MDOF  system 
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Figure  9:  Simulated  response  of  the  MDOF  system  to  the 
filtered  pseudo-random  input. 
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integration  algorithms.  The  discrete  time  transfer  function  in  the  form  of  an  ARMA  model 
is  a  numerical  integration  scheme  requiring  similar  sampling  rates.  The  parameters  of  the 
ARMA  model  are  just  weights  of  the  numerical  integration  scheme.  The  goal  of  this  report 
is  to  investigate  methods  for  determining  these  weights  from  real  experimental  structural 
response  data,  so  that  the  ARMA  model  can  be  determined  from  experimental 
measurements  and  then  used  to  predict  vibration  response.  The  following  sections  describe 
numerical  methods  that  can  be  used  to  estimate  the  parameters  of  the  ARMA  model. 
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PARAMETER  IDENTIFICATION 


The  goal  of  system  modeling  and  parameter  identification  is  to  decide  on  the  form 
of  the  model  and  then  identify  the  parameters  of  the  model.  Often  the  choice  of  the  model 
form  and  parameterization  of  the  model  are  done  at  the  same  time.  Here,  the  focus  is  on  the 
use  of  the  discrete  time  transfer  function  or  the  ARMA  model  for  structural  vibratory 
prediction.  The  model  form  is  fixed  except  for  the  choice  of  the  size  (or  order)  of  the 
model.  Experimental  data  from  the  structure  itself  are  needed  to  determine  the  proper 
order.  Once  the  order  is  selected,  experimental  data  are  used  to  parameterized  the  ARMA 
model.  Some  consideration  must  be  made  of  the  choice  of  experimental  data  and  the  data 
acquisition  procedures.  Once  data  are  collected,  an  identification  algorithm  must  be  chosen 
to  parameterize  the  model.  This  section  describes  some  of  the  difficulties  in  data 
acquisition  procedures  and  model  order  selection  based  upon  the  free  response 
autocorrelation  matrix. 

Data  Acquisition  and  Filtering 

In  order  to  parameterize  a  SISO  ARMA  model,  one  needs  to  simultaneously 
measure  the  input  and  output  time  series.  Usually  the  data  is  acquired  with  the  aid  of  anti¬ 
aliasing  filters.  The  anti-aliasing  filters  are  low  pass  filters  which  are  used  to  limit  the 
bandwidth  of  the  experimental  signals.  They  are  helpful  in  eliminating  high  frequency  noise 
and  high  frequency  modes.  The  filters  truncate  a  high  or  infinite  degree  of  freedom  system 
to  a  lower  order  system  with  modal  frequencies  in  the  allowable,  selectable  bandwidth.  An 
ideal  low  pass  filter  causes  no  phase  or  magnitude  distortion  of  the  sinusoidal  components 
of  the  signal  below  the  cutoff  frequency  and  truncates  the  sinusoidal  components  of 
frequencies  above  the  cutoff  frequency.  Practical  low  pass  filters  distort  the  phase  and 
magnitudes  of  the  sinusoidal  components  of  the  signal  below  the  cutoff  frequency  and  only 
attenuate  the  higher  frequency  components.  Lx)w  pass  filters  distort  the  signal  frequency 
according  to  the  filter's  transfer  function.  The  Fourier  transform  of  the  signal  prior  to 
filtering  (Xin(co))  is  altered  to  yield, 

^out  ~  F(w)  ^in  (24) 

where  F(o))  is  the  frequency  response  of  the  low  pass  filter.  The  z-transfonn  is 
analogously  altered  by  the  filter 


•^out  (z)  =  F(z)  Xjn  (7.) 


(25) 
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where  F(z)  is  the  discrete  time  equi’’.alent  of  F(q)).  When  both  the  measured  input  and  the 
output  time  series  have  been  conditioned  by  filters,  the  discrete  time  transfer  function 
between  the  measured  input  and  the  output  time  function  is 

Y(z)  F  (z) 

G(z)  = - (26) 

U(z)  F^iz) 

where  Fu(z)  and  Fy(z)  are  the  discrete  time  transfer  functions  of  the  input  and  output  anti¬ 
aliasing  filters.  The  discrete  time  transfer  function,  G(z),  will  be  unaffected  by  the  filters 
when  the  two  anti-aliasing  filters  have  the  same  frequency  response  function.  Thus,  if 
filters  are  identical  and  are  set  at  the  same  cutoff  frequency,  the  identification  procedure  can 
treat  the  measured  response  as  if  it  had  not  been  filtered. 

If  the  filters  are  not  identical,  then  the  measured  data  must  be  corrected  so  that  the 
filters  do  not  affect  the  transfer  function.  The  correction  takes  place  in  the  frequency 
domain  using  the  convolution  properties  of  the  Discrete  Fourier  Transform  (DFT) 
(reference  12).  The  measured  time  series  is  first  converted  into  the  frequency  domain  using 
an  FFT.  The  DFT  of  the  signal  is  corrected  to  compensate  for  the  filter  and  the  corrected 
frequency  domain  response  is  converted  back  into  the  time  domain  using  an  inverse-FFT. 
The  difficulties  with  the  procedure  is  that  the  DFTs  have  leakage  and  aliasing  errors  that 
may  result  in  a  further  noise  corruption  of  measured  data.  Also  the  transfer  function  of  the 
filter  must  be  known  precisely  to  compensate  in  the  frequency  domain. 

There  is  also  resolution  noise  introduced  by  the  digitization  of  the  signal.  The 
analog-to-digital  converters  have  a  finite  number  of  values  to  assign  the  infinite  possibilities 
of  a  voltage  range.  For  example,  4096  values  may  be  assigned  by  the  computer  to  voltages 
between  +1  and  -1  Volt.  If  only  10  percent  of  the  voltage  range  is  used,  only  410 
different  values  would  be  used  to  define  the  continuously  varying  signal.  The  resolution 
noise  is  reduced  by  using  as  much  of  the  entire  voltage  range  as  possible,  thus  amplifiers 
may  be  needed.  The  amplifiers  may  introduce  more  noise  but,  one  hopes,  not  as  much  as 
the  resolution  noise.  Also,  the  amplifiers  may  introduce  additional  distortion  of  the  signal, 
since  the  amplifiers  have  a  transfer  function  of  their  own.  The  goal  of  the  data  acquisition 
process  is  to  digitize  the  input  and  output  time  functions  while  minimizing  the  noise  and 
distortion. 
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System  Order 

Once  data  are  taken,  the  ARMA  model  can  be  parameterized.  The  order  of  the 
ARMA  model  has  to  be  selected  for  a  given  system.  The  order  of  the  ARMA  model  given 
in  Eqn  14  is  the  size  of  the  integer  p.  If  the  system  order  is  unknown  apriori,  it  has  to  be 
arbirarily  selected  or  determined  from  experimental  data.  The  free  response  can  be  used  to 
determine  order.  Consider  what  happens  to  Eqn  14  when  there  is  no  input:  the  model 
reduces  to 


y(k)  =  -ajy(k-l)  -  a2y(k-2)  - ...  -apy(k-p)  (27) 

where  the  parameters,  the  a's,  are  known  as  the  Autoregressive  parameters  of  the  ARMA 
model.  Equation  27  can  be  put  in  the  form 

[zP  +  ajzP‘4  ...+  ap  ,z  +  apl  y(k)  =  C(z)  y(k)  =0  (28) 

Processes,  like  y(k),  that  satisfy  Eqn  22  are  called  Autoregressive  (AR)  processes.  The 
model  given  by  Eqn  27  is  called  an  AR  model.  Note  that  the  order  of  the  AR  model  for 
system  free  response  is  identical  to  the  order  of  the  ARMA  model  for  system  forced 
response. 

A  summation  of  damped  sinusoids,  the  free  response  of  linear  structural  systems, 
can  be  modeled  exactly  as  an  AR  process.  (See  Appendix  B.)  Furthermore,  the  order  of 
the  AR  model  is  twice  the  number  of  modal  frequencies  or  number  of  degrees  of  freedom 
in  the  system.  The  modal  frequencies  and  damping  factors  can  be  extracted  from  the  roots 
of  C(z),  termed  the  characteristic  equation  in  the  time  domain. 

A  crude  procedure  for  determining  system  order  is  to  count  the  number  of  peaks  or 
modes  in  the  the  power  spectrum  of  the  system’s  free  response  or  count  the  number  of 
peaks  in  the  system's  frequency  response  plot  (often  automatically  displaced  on  a  spectrum 
analyzer).  Problems  occur  when  the  signal  is  too  noisy,  the  modes  are  too  closely  space  to 
distinguish  separate  peaks,  or  one  or  more  of  the  modes  is  highly  damped,  reducing  the 
size  of  one  or  more  peaks.  Though  the  procedure  can  be  effective  for  selecting  order,  it  is 
not  well-suited  for  digital  automation.  A  method  better  suited  for  automation  depends 
upon  the  discrete  time  free  response  autocorrelation  matrix.  Consider  an  AR  process  of 
order  p  given  by  Eqn  28  rearranged  as 
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(29) 


y(m-k)  =0  ,  =  1. 


If  Eqn  29  is  multiplied  by  y(q-m)  and  the  expected  value  taken,  the  result  is  the  well  known 
autocorrelation  relation  (reference  1 3) 


a^Ryy(q+k)  =0 


The  autocorrelation  matrix,  arbitrarily  dimensioned  mxm. 


Ryy(0)...Ryy(-m) 
Ryy(m)...  Ryy(O) 


(30) 


(31) 


only  has  p  linearly  independent  columns  (if  m>p).  Thus  it  is  possible  to  determine  system 
order  by  determining  the  rank  of  the  autocorrelation  matrix.  The  only  difficulty  is  in 
accurately  estimating  the  autocorrelation  matrix.  Consider  the  measured  time  series 

x(k)  =  y(k)  +n(k)  (32) 

where  y  is  the  signal,  n  is  the  noise,  and  x  is  the  measured  signal.  The  autocorrelation 
matrix  of  the  signal  is  estimated  by  the  autocorrelation  matrix  of  the  measured  signal  whose 
autocorrelation  matrix  is 

1  N-m 

^^xlj  (33) 


where  N  is  the  number  of  data  points  and  m  is  the  size  of  the  autocorrelation  matrix.  The 
rank  of  a  matrix  can  be  evaluated  using  its  singular  value  decomposition.  (See  Appendix 
C.)  A  matrix  of  size  m  and  rank  p  will  have  m-p  zero  singular  values.  The  autocorrelation 
matrix  of  y(k)  will  have  rank  p,  but  the  matrix  is  unavailable  since  y(k)  is  not  a  measurable 
value.  The  autocorrelation  of  matrix  of  the  measured  signal  x(k)  is  available  and  will  have 
full  rank  due  to  effect  of  the  noise  present  in  the  signal.  Consider  the  expected  value  of 
Rxx  when  the  noi.se  process  is  statistically  independent  from  the  signal  y(k)  and  is  an 
independent  identically  distributed  process  (reference  14)  with  variance 
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E  [R.J  =  R„  +  0^1 


(34) 


An  independent  identically  distributed  process,  n(t),  has  statistically  independent 
realizations  at  different  times  but  each  realization  comes  from  the  same  parent  distribution. 
The  rank  of  Rxx  is  thus  full  and  not  equal  to  the  rank  of  Ryy.  Some  properties  of  the 

Singular  Value  Decomposition  (SVD)  and  symmetric  matrices  can  be  used  to  reduce  Eqn 
34  to  a  useable  form.  The  Singular  Value  Decomposition  of  Ryy^  since  it  is  symmetric,  is 

Ryy=UZ  UT  (35) 

where  U  is  an  orthonormal  matrix.  The  identity  matrix  can  be  represented  by  UU^, 
reducing  Eqn  34  to 

E[R,J  =  U(I:  +  o2i)u'^  (36) 

The  rank  of  Ryy  can  be  determined  by  counting  the  number  of  singular  values  of  Rxx 
above  the  variance  of  the  noise.  There  are  at  least  three  possible  problems  with  such  a 
procedure.  First,  the  variance  of  the  noise  is  unknown.  Second,  Eqn  36  represents  the 
expected  value  of  Rxx  which  can  only  be  estimated.  Third,  the  noise  may  not  be  an 

independent  identically  distributed  process.  In  fact,  the  expected  value  of  the  noise 
autocorrelation  may  not  be  a  diagonal  matrix.  Though  there  are  some  limitations  in  the  use 
of  Eqn  36  to  determine  the  system's  order,  it  is  a  useful  relation  which  has  proven  helpful 
in  evaluating  system  order. 

Examples 

The  expected  value  of  the  autocorrelation  matrix  for  a  independent  identically 
distributed  noise  process  is  equal  to  the  variance  times  the  identity  matrix.  One  expects  a 
flat  singular  value  "spectrum"  of  the  autocorrelation  matrix.  The  singular  value  (SV) 
spectrum  is  merely  a  plot  of  the  singular  values  of  a  matrix  (in  decreasing  order)  against  the 
index  of  the  singular  values.  Figure  10  shows  the  singular  value  spectra  of  unity  variance 
Gaussian  noise  for  a  varying  number  of  points  developed  using  Eqn  33.  The  size  of  the 
autocorrelation  matrix  (20x20)  was  arbitrarily  chosen.  The  figure  shows  that  the  spectrum 
becomes  flatter  as  the  number  of  points  used  to  estimate  the  autocorrelation  matrix 
increases.  When  the  noise  SV  spectrum  is  not  flat,  it  is  roughly  linear  with  the  central 
signal  value  (in  this  case,  value  number  10)  roughly  equal  to  the  variance  of  the  noise. 
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Figure  10;  Singular  Value  Spectra  of  the  autocorrelation  matrix 
of  unity  variance  gaussian  noise. 


32 


Figure  1 1  shows  the  spectrum  for  a  noiseless  6th  order  AR  process  (512  pts).  The  signal 
is  composed  of  three  damped  sinusoids;  it  is  therefore  a  simulated  free  response.  There 
should  be  only  6  non-zero  singular  values.  The  spectrum  does  drop  from  roughly  10  for 
the  sixth  singular  value  to  2E-14  for  the  seventh  singular  value.  The  seventh  to  the 
twentieth  singular  values  range  from  2E-14  to  2E-15,  approximately  equal  to  zero  on  the 
computer  used  for  these  calculations.  Figure  11  indicates  that  the  rank  of  the 
autocorrelation  matrix  is  six. 

Determining  the  rank  of  the  autocorrelation  matrix  when  there  is  no  noise  is 
relatively  simple,  but  when  noise  is  present  in  the  signal,  Eqn  36  can  be  used.  Equation  36 
says  that  the  singular  value  spectrum  of  a  noisy  signal  is  the  addition  of  the  noise  and  the 
signal  spectra.  Figure  12  shows  the  spectra  of  the  the  noise,  the  signal,  and  the  noise  plus 
the  signal.  The  figure  demonstrates  that  Eqn  36  is  approximately  true.  With  the 
knowledge  that  the  noise  spectrum  is  roughly  linear,  one  would  be  able  to  surmise  that  the 
order  of  the  system  is  six.  Compare  the  ease  of  order  determination  using  the  singular 
value  spectrum  versus  the  use  of  the  Power  Spectrum  shown  in  Figure  13.  A  structural 
dynamicist  would  be  able  to  recognize  that  there  are  three  significant  peaks,  but  automating 
such  a  decision  might  be  rather  difficult.  The  problem  is  aggravated  when  there  are  closely 
spaced  modes.  Figure  14  shows  the  singular  value  spectra  of  a  sixth  order  system  with  a 
pair  of  closely  spaced  modes  (specifically  at  14.6  and  15.0  Hz,  which  are  roughly  one  DFT 
frequency  bin  apart)  when  the  random  noise  is  added.  Extending  the  linear  noise  portion  of 
the  spectrum  back  indicates  that  a  sixth  order  model  is  appropriate.  The  power  spectrum 
shown  in  Figure  15  is  ambiguous.  It  is  very  difficult  to  tell  if  there  are  three  modes  present 
or  if  the  spectrum  is  noisy.  The  singular  value  spectrum  may  not  be  a  "foolproof  method 
to  select  the  system  order,  but  the  decision  process  when  using  the  singular  value  spectrum 
is  less  subjective  and  more  automated.  If  there  is  some  doubt,  the  higher  order  model 
should  be  used  because  overspecification  of  the  system  order  is  not  disastrous.  An 
overspecified  model  may  perform  as  well  as  a  model  with  the  proper  order,  but  an 
underspecified  model  will  not. 

Overview 

The  problems  of  data  acquisition  and  order  selection  were  briefly  discussed  in  this 
section.  Free  response  can  be  used  to  determine  system  order  through  inspection  of  the 
singular  value  spectra.  It  must  be  noted  that  often  it  may  be  difficult  or  impossible  to  obtain 
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Figure  11:  Singular  Value  Spectra  of  the  autocorrelation  matrix 
of  a  signal  containing  three  damped  sinusoids. 
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Figure  15:  Power  Spectrum  of  a  noisy  signal. 
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free  response  or  the  impulse  response  data.  The  system  order  may  have  to  be  selected 
while  parameterizing  the  ARMA  model.  This  can  be  accomplished  by  increasing  the  order 
of  the  ARMA  model  during  the  identification  procedure  until  the  estimated  model  is 
acceptable,  by  some  measure.  The  resulting  order  of  the  ARMA  model  may  not  be  the 
same  as  it  would  be  using  free  response,  and  the  ARMA  model  may  end  up  being 
excessively  large.  Even  if  the  free  or  impulse  response  is  available,  the  proper  order  may 
not  be  selected  as  shown  by  the  examples.  Given  an  estimate  of  the  ARMA  model  order 
and  experimental  data,  the  next  step  is  to  estimate  the  ARMA  model.  The  next  section 
presents  identification  procedures  to  do  the  ARMA  parameter  estimation. 
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IDENTIFICATION  OF  THE  ARMA  MODEL 


This  section  presents  some  classical  time  domain  identification  methods  that  use 
experimental  data  to  estimate  the  parameters  of  the  ARMA  model.  The  experimental  data  is 
used  to  write  an  overdetermined  set  of  equations.  The  most  straightforward  way  of  solving  * 

these  equations  is  the  Least  Squares  method.  Unfortunately,  the  Least  Squares  method 
produces  biased  parameter  estimates.  Many  other  methods  exist  that  attempt  to  reduce  the 
bias  in  some  fashion.  Some  of  these  methods  are  included  here. 

The  ARMA  model 

The  ARMA  model  contains  transient  response  terms,  the  AR  part  of  the  model,  and 
forced  response  terms,  the  MA  part  of  the  model.  The  model  can  be  written  in  vector  form 
as 

- -^k-p)  4^(k)  %(k-l)  4^(k-p)]0  +  e(k)  (37) 


where  the  subscript  m  denotes  measured  quantities,  e(k)  denotes  a  residual  or  error,  and 

e=(ai,  a2, ap,  bQ,  bj . bp]T  (38) 

Identification  algorithms  are  used  to  estimate  the  parameter  vector,  0,  by  some  criteria 
given  multiple  equations  produced  by  varying  the  index  k  in  Eqn  37.  In  order  to  identify 
the  parameters  of  the  ARMA  model,  Eqn  37  must  be  written  for  at  least  2p+ 1  different  time 
instances.  If  the  data  is  exact  and  from  an  ARMA  process,  the  exact  parameters  can  be 
found  by  solving  the  2p+l  equations.  When  noise  corrupts  the  measured  input  and  output 
discrete  time  histories,  it  is  impossible  to  determine  the  ARMA  model  parameters  exactly. 
An  estimate  of  the  parameters  can  be  made  by  solving  an  overdetermined  set  of  equations 
(more  than  2p+l  equations).  There  are  various  methods  to  solve  this  set  of  equations.  The 
methods  have  varying  criteria  and  will  be  presented  in  the  following  sections. 

The  Least  Squares  (LSI  algorithm 

The  Least  Squares  algorithm  solves  the  overdetermined  set  of  equations  by 
minimizing  the  sum  of  the  squares  of  the  residuals.  If  there  are  N  data  points  in  each  time 
history,  N-p  equations  can  be  written,  where  p  is  the  model  order.  The  equations  in  matrix 
form  are 
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“8m  ® 


(39) 


where 

8m=[ym(P+l).  ym(P+2) . ym(N)]'^ 

(40) 

* 

and 

e=[e(p-H),  e(p+2)....,  e(N)]T 

(41) 

-ym(P)  -  -ym^l) 

m 

-ym(p+l)-  -ym(2)  u^(p-^2)...  u^(2) 

(42) 

_-y^(N-l)...-y^(N-p)  u^(N)  ...u^(N-p)_ 

Solving  Eqn  39  by  minimizing  the  sum  of  the  squares  of  the  residual  is  equivalent  to 
solving  a  set  of  equations  termed  the  normal  equations, 

(Am^Am)  6  LS  ^m^8m  (“^3) 

The  solution  of  the  normal  equations  is 

(^m^^m)*^  ^m^Bm  (‘^) 

where  the  term  (Am'^Am)*^  is  referred  to  as  the  pseudo-inverse  of  Aj^.  The  normal 
equations  become  highly  ill-conditioned  as  either  the  sample  rate  or  the  number  of 
parameters  increase.  One  method  to  reduce  the  ill-conditioning  problem  is  to  solve  the 
overdetermined  set  of  equations  through  a  technique  termed  singular  value  decomposition. 
Solving  the  equations  by  singular  value  decomposition  will  result  in  a  LS  estimate,  but  the 
equations  are  solved  in  a  different  fashion.  (See  Appendix  D.) 

The  bias  in  the  Least  Squares  estimate  is  apparent  when  Eqn  39  is  premultiplied  by 
the  the  pseudo-inverse  of  Aj^,  resulting  in 


®LS  -  9  -  ^  A^^'g 


(45) 


Taking  the  expectation  of  Eqn  45,  yields 
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E  [BlS  ]  =  0  -  E  [(A^TAm)- 1  A^Te]  (46) 

Equation  45  shows  that  the  LS  parameter  estimates  are  biased  by  the  E[(ArjiTAjyj)-l 

Am^e]  term.  It  has  been  widely  recognized  that  the  LS  estimates  are  biased  (references  15 

and  16).  This  sensitivity  of  the  LS  estimates  to  noise  often  makes  the  estimates  unreliable.  * 

Usually  an  overspecification  of  the  model  order,  which  will  be  discussed  later,  or  an 

iterative  technique  to  improve  the  estimates  can  be  used  to  reduce  the  bias.  » 

Generalized  Least  Squares 

The  Generalized  Least  Squares  (GLS)  algorithm  is  an  iterative  technique  used  to 
improve  the  biased  LS  parameter  estimates.  The  algorithm  attempts  to  reduce  the  noise 
bias  by  modelling  the  residuals.  The  residual  process  is  modeled  as  an  Autoregressive 
(AR)  process 

0  =  C(z)  e(k)  -  e(k)  (47) 

where  e(k)  is  another  residual  process.  The  model  of  the  system  is  an  ARMA  model 

A(z)  ym(k)  =  B(z)  u(k)  +  e(k).  (48) 

Multiplying  Eqn  48  by  C(z)  and  using  Eqn  47,  yields 

A(z)  C(z)  ym(k)  =  B(z)  C(z)  Um(k)  +  E(k).  (49) 

If  one  defines  C(z)  ym  (k)  as  x(k)  and  C(z)  Um(k)  as  z(k),  then  the  new  system  model  is 

A(z)  x(k)  =  B(z)  z(k)  +e(k)  (50) 

which  is  the  basis  for  the  GLS  algorithm.  The  GLS  estimate  of  the  system  model  is  the  LS 
estimate  using  the  transformed  processes  x(k)  and  z(k).  The  GLS  algorithm  in  step  form 
is;  • 

1.  Estimate  9ls  (  Eqn  44), 

2.  Calculate  the  residual  sequence  e(k). 

3.  Model  the  residual  sequence  as  an  Autoregressive  Process. 

4.  Generate  x(k)  and  z(k). 

5.  Compute  0GLS-  ES  estimate  of  (Eqn  50) . 

6.  Go  to  Step  (2)  until  OgLS  converges. 
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The  GLS  algorithm  replaces  the  residual  process  of  the  LS  estimate,  e(k),  by  the  residual 
process  of  the  Generalized  Least  Squares  estimate.  The  bias  of  the  GLS  estimate  is 
hopefully  smaller  than  the  bias  of  the  LS  estimate  because  e(k)  will  be  of  smaller 
magnitude  than  e(k). 

Reference  17  proposed  a  modification  of  the  GLS  algorithm  called  Modified  Least 
Squares  (MLS)  estimation  algorithm.  The  MLS  algorithm  also  models  the  residual  process 
as  an  Autoregressive  process.  The  algorithm  is  noteworthy  because  it  readily  demonstrates 
the  effect  of  the  residual  model.  Equation  44  is  rewritten  as 

0  =  0LS  +  (51) 

If  e  is  estimated  as  an  AR  process,  then  there  is  a  set  of  overdetermined  equations  based 
upon  the  residual  sequence 

e  =  Q  c  +  e  (52) 

The  residual  of  Eqn  52,  e,  is  neglected  and  the  Modified  Least  Squares  estimate  is  written 
from  Eqn  51  as 

0  MLS=  0LS  +  (Am^Am)-  ^  A^T  q  q  (53) 

The  Generalized  Least  Squares  algorithms  model  the  bias  term  and  add  the  result  to  the  LS 

i 

estimate  in  order  to  improve  the  estimate.  The  difficulties  of  the  GLS  algorithms  are  that 
they  are  iterative  (and  thus  prone  to  computational  stability  concerns)  and  the  size  of  the 
residual  model  has  to  be  specified. 

Instrumental  Variables 

The  GLS  method  relies  upon  the  AR  model  of  the  residuals  to,  in  a  sense,  estimate 
the  bias  of  the  LS  estimates  and  therefore  better  estimate  the  parameters.  The  Instrumental 
Variable  (IV)  technique  attempts  to  force  the  bias  towards  zero  by  proper  selection  of  a 
different  "pseudo-inverse."  Suppose  that  a  matrix,  Z,  could  be  found  so  that 

E[ZTe]  =0 

and  E[  Z^Aml  =  R  (54) 

where  R  is  a  non-singular  matrix.  The  "pseudo-inverse"  of  (Z^Afjj  can  be  formed. 
Premultiplying  Eqn  41  by  this  "pseudo-inverse"  yields 
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eiv  =  0-(ZTAm)-l  ZTe 


(55) 


where  Biv  =  (Z'^Am)'^  Z^gj^  (56) 

Here  the  expected  value  of  the  IV  estimate  is  the  parameter  vector.  The  main  problem  with 
the  IV  method  is  the  selection  of  Z  .  The  matrix  Z  is  comprised  of  instrumental  variables 
and  is  structured  similar  to  A^.  The  instrumental  variables  should  be  correlated  to  the 

noise-free  signal  and  uncorrelated  to  the  noise.  The  perfect  set  of  instrumental  variables  is 
the  noise-free  signal  itself  which  is  unavailable.  An  auxiliary  model  can  be  used  to  generate 
the  instrumental  variables  and  hence  the  matrix  Z.  Reference  18  has  shown  that  any  stable 
model  of  the  same  order  as  the  system  provides  consistent  unbiased  IV  estimates,  when 
there  is  no  noise  in  the  input  measurements.  The  Z  matrix  is  constructed  from  the 
instrumental  variables  by  substituting  Z  for  Aj^j  and  the  output  of  the  auxiliary  model  for 

the  output  of  the  system.  Reference  15  has  also  assumed  that  the  input  is  free  of  noise.  The 
assumption  that  the  input  is  noise  free  may  be  valid  in  control  engineering  where  the  input 
can  be  a  reference  voltage  but  the  input  will  surely  contain  noise  in  structural  applications. 

Maximum  Likelihood  Method 

The  Maximum  Likelihood  Method  (MLE),  like  the  GLS  algorithm,  assumes  that  the 
residual  sequence  is  an  AR  process  and  uses  the  model 

bjUn,(k-i)  -t-  ^Cje(k-i)  (57) 

i=o 

The  parameters,  the  aj,  bj,  and  Cj,  are  estimated  to  maximize  the  probability  of  the  e(k) 
sequence  under  the  normal  probability  density  function.  References  3  and  19  give  an 
iterative  algorithm  to  perform  the  estimation.  The  major  difficulties  with  the  MLE  method 
are  that  the  method  is  iterative  and  the  assumption  of  normal  residuals  may  not  be  valid  in 
many  cases. 


m  =  ]^a.yjk-i)  +  ^ 

i=0 


i=0 


The  Overspecified  LS  model 

The  main  difficulty  with  the  time  domain  estimation  of  the  ARMA  model  is  the 
estimation  bias.  The  LS  algorithm  has  been  modified  to  produce  the  other  major  algorithms 
in  hopes  of  reducing  or  eliminating  the  bias.  The  major  difficulties  of  the  algorithms  are 


44 


that  they  are  iterative  in  nature  and  lack  the  simplicity  of  the  LS  algorithm.  In  order  to 
implement  the  LS  algorithm  one  needs  to  "measure"  the  time  series  data,  form  matrices  of 
the  proper  order,  and  invert  a  single  matrix.  There  are  no  convergence  considerations.  An 
alternative  way  to  reduce  the  bias  is  to  overspecify  the  LS  model.  The  additional 
parameters  serve  as  additional  modes  which  may  be  used  to  effectively  model  the  noise.  If 
the  noise  is  adequately  modeled,  the  residuals  are  reduced  and  the  estimation  bias  is 
reduced.  However,  the  overspecified  model  will  contain  extra  poles  and  zeroes  related  to 
the  noise.  If  the  overspecified  model  is  used  to  predict  vibration,  the  noise  modes  will  still 
be  present.  Usually  the  noise  modes  in  the  overspecified  model  are  of  such  a  low 
magnitude  that  the  prediction  is  unaffected. 

Examples 

Four  computer  programs  based  on  the  LS,  the  GLS  (actually  the  MLS  form  of  the 
GLS  algorithm),  the  IV  (as  proposed  in  reference  18),  and  the  MLE  (according  to  reference 
19)  were  completed  and  simulated  data  contaminated  with  noise  was  used  to  test  the 
algorithms.  The  Laplace  domain  model  of  the  test  system  was 

G(s)=(2s+100)  /  s2+2s+100)  (58) 

which  was  transformed  by  Tustin's  Bilinear  rule  (sampled  at  4(X)  rad/sec)  to 

0.02141  z^+  0.01207  z  -  0.009335 

G(z)  - - » -  (59) 

z  -  1.945  z  +  0.969 


in  the  discrete  time  domain.  Input  time  histories  were  developed  to  nxxlel  random  input  to 
the  system.  Random  numbers  (uniform  distribution)  were  determined  to  represent  the  input 
time  increments.  The  magnitude  of  the  input  was  gaussian.  If  the  input  time  increment  was 
greater  that  the  system  sample  time,  the  input  was  assumed  to  vary  linearly  between  input 
specification  times.  This  was  a  simple  attempt  to  set  the  "wave"  length  and  spectrum  of  the 
input  time  series.  The  random  time  series  was  used  as  input  to  Eqn  59  to  produce  a  noise- 
free  system  output  time  series.  The  output  was  then  contaminated  with  simulated  Gaussian 
noise  at  an  signal-to-noise  ratio  (SNR)  of  40  dB.  The  signal-to-noise  ratio  is  defined  as 


SNR  =  20  logjQ 


\ms 

^noise 


(60) 
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where  is  the  root  mean  square  amplitude  of  the  noise-free  signal.  Figures  16,  17, 
and  18  show  three  different  random  inputs  that  were  used  in  the  numerical  experiments. 
Table  1  shows  the  parameter  estimates  of  the  four  algorithms  for  each  case. 

It  is  clear  that  the  IV  method  gives  the  best  estimates  of  the  AR  parameters.  The 
MA  parameters  are  difficult  to  obtain  and  the  estimation  success  is  dependent  upon  the 
excitation  (the  input).  The  third  input,  which  has  the  widest  excitation  range,  appears  to 
provide  the  best  parameter  estimates.  This  is  related  to  the  concept  of  persistent  excitation 
and  has  not  been  discussed  as  yet.  A  persistent  excitation  is  an  input  that  provides  enough 
information  so  that  the  Moving  Average  parameters  can  be  accurately  estimated.  Though 
all  three  inputs  appear  to  be  persistently  exciting,  there  seems  to  be  an  optimal  input  among 
the  three.  The  concept  of  optimal  input  is  discussed  in  reference  20. 

The  LS  estimates  are  biased  when  noise  is  present.  System  parameters  developed 
using  a  non-overspecified  LS  approach  would  not  be  suitable  for  predicting  vibratory 
response.  Among  the  three  iterative  methods,  the  IV  method  was  the  most  accurate 
method,  the  easiest  to  program,  the  quickest  computationally,  and  the  most  stable  method. 
The  main  difficulty  with  the  IV  algorithm  is  the  generation  of  the  instrumental  variable. 
The  method  suggested  in  reference  15  requires  that  the  input  must  be  noise-free  to  insure 
convergence.  If  the  input  is  a  control  system  reference  voltage,  the  assumption  of  noise- 
free  input  is  valid.  If  the  input  is  a  structural  measurement,  the  assumption  is  invalid. 

Overview 

A  number  of  parameter  identification  algorithms  have  been  discussed.  All  the 
algorithms  are  schemes  to  solve  an  overdetermined  set  of  equations  in  order  to  estimate  the 
parameters  of  a  model.  The  key  difficulty  to  be  overcome  is  the  effect  of  noise  bias.  Many 
of  the  algorithms  utilize  an  iterative  scheme  to  model  the  noise  and  improve  the  parameters. 
Other  non-iterative  algorithms  simply  overspecify  the  model  order  and  use  the  additional 
modes  of  the  model  as  a  noise  outlet.  Usually  the  overspecified  algorithms  out-perform  the 
iterative  schemes  both  in  computer  cost  and  accuracy.  The  difficulty  with  the 
overspecified  algorithms  is  the  sorting  of  the  noise  modes  from  the  true  modes  of  the 
model.  Time  domain  modal  techniques  such  as  the  ITD  and  the  ERA  also  use 
overspecification  to  reduce  noise  bias.  These  techniques  use  modal  confidence  factors 
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Input  (9) 


Figure  16:  Input  Case  1. 
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Figure  17:  Input  Case  2. 
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Figure  18:  Input  Case  3. 
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Table  1:  Parameter  Estimates 


al 

a2 

bO 

bl 

b2 

True  values: 

-1.945 

0.969 

0.0214 

0.0121 

-0.00934 

Input  =  Case  1 

LS 

-1,740 

0.769 

0.148 

-0.116 

0.00016 

GLS 

-1.875 

0.892 

0.135 

-0.143 

0.0238 

rv 

-1.947 

0.971 

0.128 

-0.196 

0.0929 

MLE 

-2.802 

1.654 

0.383 

-0.712 

0.176 

Input  =  Case  2 

LS 

-1.829 

0.858 

0.040 

-0.024 

0.0190 

GLS 

-1.860 

0.884 

0.034 

-0.012 

0.0085 

IV 

-1.942 

0.967 

0.040 

-0.022 

0.0094 

MLE 

-2.674 

1.745 

0.064 

-0.070 

-0.0307 

Input=  Case  3 

LS 

-1.788 

0.812 

0.0194 

0.0154 

-0.0041 

GLS 

-1.857 

0.877 

0.0208 

0.0138 

-0.0059 

rv 

-1.941 

0.966 

0.0197 

0.0128 

-0.0082 

MLE 

-1.924 

0.949 

0.0207 

0.0116 

-0.0087 
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developed  from  mode  shape  information  to  sort  the  noise  modes.  The  simplest  and  most 
efficient  technique  to  identify  an  SISO  ARMA  model  is  to  overspecify  the  model  order 
before  using  the  LS  algorithm.  But  since  there  is  no  method  at  present  to  reduce  the 
overspecified  model  once  it  has  been  estimated,  the  predictive  model  would  be  excessively 
large  and  contain  noise-related  modes. 
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FREE  RESPONSE  METHODS 


The  ARMA  model  can  be  identified  in  a  single  stage  using  the  algorithms  discussed 
in  the  previous  sections  or  in  two  separate  stages.  In  a  two-stage  method,  the  AR 
parameters  are  identified  separately  from  the  MA  parameters.  The  AR  parameters  can  be 
identified  from  free  or  impulse  response  records.  Also,  if  the  the  frequencies  and  damping 
factors  are  known,  the  AR  parameters  can  be  indirectly  identified.  All  of  the  algorithms  in 
the  preceding  section  can  be  modified  slightly  to  identify  the  AR  model  from  free  response. 
In  addition,  there  are  many  time  domain  procedures  that  are  used  to  estimate  structural 
frequencies  and  damping  factors.  Among  these  arc  the  Ibrahim  Time  Domain  Technique 
(ITD),  the  Eigensystem  Realization  Algorithm  (ERA),  and  the  Complex  Exponential 
Method  (CEM).  Essentially,  these  algorithms  identify  the  pole  locations  in  the  z-plane. 
The  poles  locations  can  be  used  to  identify  an  AR  model.  The  algorithms  can  also  be  used 
to  identify  the  complex  mode  shapes  of  the  structure.  These  methods  require  free  response 
or  impulse  response  of  the  structure  measured  at  many  locations  on  the  structure.  Ideally, 
the  AR  model  can  be  identified  from  a  single  signal.  The  structural  identification 
techniques  require  more  data  and  produce  more  information  than  might  be  desirable  for  the 
transient  response  prediction  problem  which  this  report  addresses.  There  exists  many 
algorithms  in  the  signal  processing  field  such  as  Prony's  method,  the  Pisarenko  Method, 
and  the  Principal  Eigenvectors  Method  (PEM)  which  identify  the  frequencies  and  damping 
factors  from  a  single  response  record.  Recently,  these  signal  processing  tech  frques  have 
made  inroads  into  the  structural  identification  field.  The  signal  processing  and  the  lime 
domain  modal  methods  are  ail  linear  prediction  methods  (reference  13)  and  are  non¬ 
iterative.  The  major  means  of  reducing  bias  errors  using  linear  prediction  is  to  overspecify 
the  system  order.  The  overspecification  produces  extraneous  poles  which  must  be 
distinguished  from  the  signal  poles.  The  time  domain  modal  methods  are  able  to 
distinguish  the  poles  by  building  modal  confidence  factors.  The  modal  confidence  factors 
are  determined  by  comparing  the  mode  shapes  of  the  two  separate  identifications.  The 
signal  processing  techniques  do  not  identify  mode  shapes  and  the  extraneous  poles  must  be 
distinguished  from  the  signal  poles  by  some  apriori  knowledge.  The  PEM  can  be  modified 
slightly  to  automatically  determine  which  poles  correspond  to  the  signal  and  which 
correspond  to  the  noise. 

The  AR  model 

The  Autoregressive  model  given  by  Eqn  27  can  exactly  model  a  signal  comprised  of 
damped  sinusoids  (Appendix  B).  The  roots  of  the  time  domain  characteristic,  Eqn  28,  are 
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directly  related  to  the  frequencies  and  damping  factors  by  Eqn  8,  where  z  and  s  are  complex 
numbers.  The  frequencies  and  damping  factors  are  related  to  the  roots  of  the  characteristic 
equation  by 


s.=-;co±ja)j 


(61) 


when  substituted  into  Eqn  8  yields 


z.= 

1 


(cos  (O^h  +  j  sin  h) 


(62) 


Thus,  if  one  is  able  to  model  a  signal  as  an  AR  process,  a  characteristic  equation  can  be 
written  using  the  AR  parameters.  The  roots  of  the  characteristic  equation  define  the 
frequencies  and  damping  factors  of  the  sinusoids  present  in  the  signal.  Conversely,  if  the 
frequencies  and  damping  factors  are  known,  the  AR  model  can  be  determined  using  the 
approximations  outlined  above. 

The  locations  of  the  poles,  roots  of  the  characteristic  equation,  in  the  z-plane  are 
defined  by  the  frequencies,  damping  factors,  and  sample  rate  in  Eqn  62.  One  sees  that  the 
stability  region  (^  >  0)  in  the  z-plane  is  inside  the  unit  circle.  Lines  of  constant  damping  are 
decaying  spirals  and  lines  of  constant  damped  natural  frequencies  are  straight  lines  as 
shown  by  Figure  19.  According  to  Figure  19,  the  root  locations  get  closer  to  the  (1,0) 
point  in  the  z-plane  as  the  sample  rate  increases.  As  the  root  locations  get  closer  to  the 
(1,0)  point,  perturbations  in  the  root  locations  result  in  greater  errors  in  frequency  and 
damping  estimates.  Figures  20  to  25  show  the  error  envelopes  for  frequency  and  damping 
estimates  for  various  root  location  perturbations.  For  instance.  Figures  21  and  22  show 
that  the  identified  damped  natural  frequency  will  be  in  error  by  around  10%  as  the  root 
location  travels  to  the  (1,0)  point  (for  a  perturbation  radius  of  0.01)  and  the  damping 
estimate  will  be  in  error  over  400%.  The  figures  show  that  damping  error  is  much  more 
sensitive  to  the  exact  damping  value  than  is  the  frequency  estimate  and  the  size  of  the 
frequency  envelope  is  approximately  linear  with  respect  to  the  perturbation  size  (See 
Appendix  E.)  The  figures  are  different  from  those  calculated  in  reference  2.  The  figures 
illustrate  that  when  the  sample  rate  is  high  and  the  damping  factor  is  low  (as  with  most 
identified  structural  systems),  the  roots  of  the  characteristic  equation  must  be  identified  very 
accurately  for  reliable  estimates  of  frequency  and  damping  factors  and  hence  reliable 
vibratory  response  prediction.  This  means  that  the  identification  scheme  must  be  accurate 
in  the  parameterization  of  the  characteristic  equation. 
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Figure  20:  Envelope  of  percent  damped  frequency  error 
(radius  of  uncertainty  =0.01). 
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Figure  22:  Envelope  of  percent  damped  frequency  error 
(radius  of  uncertainty  =0.001). 
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Figure  24;  Envelope  of  percent  damped  frequency  error 
(radius  of  uncertainly  =0.0001). 
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Identified  Damping  Error  (%) 


Figure  25:  Envelope  of  percent  damping  error 
(radius  of  uncertainty  =0.0001). 
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Pisarenko  Harmonic  Decomposition 

The  first  method  discussed  is  the  Pisarenko  Harmonic  Decomposition  Method 
(PHD)  and  is  also  called  the  Total  Least  Squares  (TLS)  algorithm  (reference  21).  The  PHD 
method  is  based  upon  the  autocorrelation  of  an  AR  process.  Equation  30  is  written  p+1 
times  and  the  AR  coefficients  are  determined  through  an  eigensolution. 


Ryy(O) 


Ryy(-P)' 


Ryy(p)  •••  Ryy('J) 


=  0 


or 

Ryy  0  =  0  (63) 

But  Ryy  is  unavailable  and  R^x  tnust  be  used  instead.  Substitution  of  Eqn  34  into  Eqn  63 
yields 


(Rxx  •  I  )  e  =  0  (64) 

which  is  an  eigenvalue  problem.  The  parameter  vector  is  the  eigenvector  of  Rxx  which 
corresponds  to  the  eigenvalue  equal  to  the  variance  of  the  noise.  Since  the  variance  of  the 
noise  is  not  known  apriori,  one  of  the  eigenvalues  has  to  be  chosen  as  the  noise  eigenvalue. 
Typically,  the  signals  have  a  signal-to-noise  ratio  greater  than  1,  and  the  noise  eigenvalue 
will  be  the  eigenvalue  with  the  lowest  magnitude. 

The  Correlation  Fit 

Reference  21  suggests  another  methods  based  upon  the  autocorrelation  relationship. 
The  method  is  called  the  Correlation  Fit  (CF).  In  the  CF  method,  Eqn  30  is  written  more 
than  p  times 


■  Rx,(0)  ...Rxx(^-P)~ 

_Rxx(q-l)...I^x<^-p)_ 

(65) 
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with  Rxx  used  as  an  estimate  of  Ryy  The  Rxx(0)  term  in  the  autocorrelation  matrix  is  the 

most  corrupt  with  noise  because  its  expected  value  contains  the  variance  of  the  noise,  while 
the  expected  value  of  the  other  R^x  terms  are  the  corresponding  Ryy  terms.  The  CF 
eliminates  the  first  p  equations  in  Eqn  65  because  these  equations  contain  the  Rxx(O)  terms. 
The  remaining  q-p  equations  can  then  be  solved  in  a  Least  Squares  fashion.  Reference  21 
compares  the  CF  to  other  methods. 

Overspecification  of  the  AR  model 

Reference  21  also  discusses  some  of  the  advantages  of  overspecification  of  the  AR 
model.  The  process  y(t)  was  assumed  to  be  an  AR  model  of  order  p 

A(z)  y(k)  =  0  (66) 

where  A(z)  is  a  pth  order  polynomial  in  z.  Multiplying  Eqn  66  by  an  arbitrary  polynomial 
in  z  of  order  q  yields 

D(z)  A(z)  y(k)  =  C(z)  y(k)  =0  (67) 

or 

[1+Cjz'V  ...+  y(k)=0  (68) 

The  parameters  of  C(z)  can  determined  by  any  of  the  other  methods  described.  However, 
Ryy  is  singular  for  orders  greater  than  p,  so  there  aie  an  infinite  number  of  solutions  to 

Eqn  68.  But,  A(z)  is  a  factor  of  C(z),  and  therefore  the  AR  parameters  can  be  uncovered  in 
any  particular  C(z).  The  infinite  number  of  possibilities  for  C(z)  are  manifested  by  an 
infinite  number  of  possible  D(z). 

The  Principal  Eigenvectors  Method 

The  Principal  Eigenvectors  Method  (PEM)  of  identification  is  a  simple  modification 
of  an  overspecified  LS  solution.  The  PEM  finds  the  unique  overspecified  parameter  vector 
with  a  minimum  norm.  The  PEM  sets  up  the  LS  solution  from  an  overspecified  model 

®  ~  Bin  "*■  ®  (69) 

The  matrix  is  a  scaled  version  of  Rxx  and  as  such  its  rank  is  full  (from  the  noise 

contributions)  even  in  the  overspecified  form.  However,  the  first  p  singular  values 
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correspond  to  signal-plus-noise  and  the  q-p  remaining  singular  values  correspond  to  just 
noise.  The  matrix  can  be  split  by  Singular  Value  Decomposition  into  a  signal  subspace. 


(aJ  A„) 


signal 


=  >  a.  u.  V. 

1  1  1 

i=l 


and  into  a  noise  subspace 


T 

^ni^noise 


=  >  a.  u.  V 

I  I  I 


i=p+l 


(70) 


(71) 


where  the  singular  values  are  ordered  from  largest  magnitude  to  the  smallest.  Since,  a 
singular  or  ill-conditioned  solution  of  the  normal  equations  is  expected  and  the  factor  A(z) 
of  the  overspecified  model  C(z)  is  desired,  the  noise  subspace  of  can  be  zeroed 

to  reduce  the  effects  of  noise.  The  remaining  solution  is  a  singular  solution  with  the  useful 
property  of  being  a  minimum  norm  solution.  (The  length  of  the  the  overspecified  parameter 
vector  is  the  shortest  of  all  the  singular  solutions.)  The  remaining  noise  in  the  matrix  will 
affect  the  factor  D(z)  and  not  greatly  affect  A(z).  The  PEM  effectively  modifies  the 
autocorrelation  matrix  and  reduces  the  effects  of  the  noise.  The  penalty  paid  for  the 
modification  is  that  the  AR  model  had  to  be  overspecified  and  will  contain  signal  poles  and 
extraneous  poles.  The  only  remaining  difficulty  is  how  to  sort  the  roots  of  A(z)  from  the 
roots  of  D(z). 

The  Modified  PEM 

The  estimates  of  the  frequencies  and  damping  factors  are  improved  when  the  model 
is  overspecified,  but  it  is  often  difficult  to  distinguish  the  system  frequencies  and  dampings 
from  the  extraneous  ones  introduced  by  the  overspecified  model.  The  Modified  Principal 
Eigenvectors  Method  (MPEM)  for  damped  exponentials  (reference  22)  addresses  this 
problem.  The  MPEM  solves  a  set  of  LS  equations  using  the  PEM,  but,  because  of  the  way 
the  equations  have  been  written,  the  signal  poles  can  be  distinguished  from  the  extraneous 
poles.  The  MPEM  uses  the  mapping  in  Eqn  8  to  estimate  the  frequencies  and  damping 
factors.  The  singular  values  of  the  discrete  time  autocorrelation  matrix  are  obtained  in  the 
MPEM  and  the  order  of  the  system  can  be  identified  in  an  automated  fashion. 

The  MPEM  requires  expressing  the  AR  model  in  the  backward  direction  and  thus 
places  stable  poles  outside  the  unit  circle.  The  overspecified  model  is  used  to  increase  the 
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accuracy  of  the  estimation  in  the  presence  of  noise.  The  overdetermined  equations  for  a  Lth 
order  backward  linear  prediction  can  be  written  in  matrix  form  as 


ym(2)  ym(3)  •••ym(L+l) 

ym(3)  ym(4)  •••ym(L+2) 

^  ym(i)  ^ 
ym(2) 

+ 

r  e(l)  N 
e(2) 

_y^(N-L-l)3{„(N-L)...  y^(N)  _ 

^ym(N-L)  ^ 

.e(N-L), 

or 

Qm  c  =  *dm  +  e  (74) 

The  backward  polynomial,  the  polynomial  corresponding  to  the  overspecified  backward 
AR  model,  can  be  written  as 

C(z)=  1+  cjz  +  C2  +  ...  +  CL  (75) 

There  are  p  roots  of  Eqn  75  located  outside  the  unit  circle  representing  the  poles  associated 
with  the  system  frequencies  and  damping  factors.  There  are  L-p  extraneous  poles.  The 
normal  equations  are 


Qm^Qm  ^  “  "Qm^^m 


(76) 


Again  the  matrix  on  the  left  hand  side,  Qm^Qm^  split  into  a  signal  subspace  and  a  noise 
subspace,  and  the  noise  subspace  is  zeroed.  This  solution  of  the  normal  equations  forces 
the  extraneous  poles  inside  of  the  unit  circle,  thereby  allowing  the  system  poles  to  be 
distinguished  from  the  extraneous  ones  (reference  22).  The  solution  is 


c  = 


J  ^ 


(77) 


where  Uj^,  V|(  are  p  principal  eigenvalues  and  eigenvectors  of  the  singular  value 
decomposition  of  Qm'^Qm-  Reference  22  has  shown  that  the  MPEM  does  achieve  the 
Cramer-Rao  lower  bound  for  the  variance  of  the  frequency  and  damping  estimates.  Once 
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the  system  frequencies  and  dampings  are  estimated  from  the  MPEM,  Eqn  8  can  be  used  to 
estimate  the  AR  parameters  at  any  sampling  rate. 

Examples 

1.  Simulated  Free  Response 

Simulated  free  response  records  were  generated  and  corrupted  with  various  kinds 
of  noise.  The  major  signal  processing  algorithms  were  used  to  estimate  the  damped  natural 
frequencies  and  damping  factors  from  the  corrupted  response.  The  types  of  corruption 
ranged  from  gaussian  noise  to  correlated  noise  to  constant  bias.  The  simulated  signals 
contained  one,  two,  or  three  sinusoidal  components  and  256  points  sampled  at  2  kHz.  The 
SDOF  simulated  data's  natural  frequency  was  1  Hz  with  1  percent  damping.  The  two  DOF 
simulated  data  contained  frequency  components  at  1  and  10  Hz,  with  1  percent  damping. 
The  three  DOF  simulated  data  contained  frequency  components  at  1,10,  and  25  Hz,  all 
damped  at  1  percent.  The  sample  rate  was  two  thousand  times  greater  than  the  first  system 
frequency  which  placed  the  pole  location  very  close  to  the  (1,0)  point.  According  to  the 
error  envelopes  discussion  at  the  beginning  of  this  section,  at  such  a  high  sampling  rate  the 
identification  routines  would  have  to  estimate  the  pole  location  very  accurately  in  order  to 
accurately  estimate  the  frequency  and  damping  factor.  In  order  to  represent  a  system  with 
natural  frequencies  such  as  these,  one  would  need  a  sample  rate  of  at  least  500  Hz  (twenty 
times  the  highest  frequency).  The  sample  rate  that  was  used  is  four  times  greater  than  500 
Hz  and  represents  an  extremely  difficult  identification  problem. 

The  algorithms  used  in  the  simulation  were  the  LS,  IV,  PHD,  MLE,  GLS,  CF, 
and  MPEM.  The  LS  algorithm  was  used  in  the  classical  form  and  in  the  overspecified 
form.  The  overspecified  LS  estimate  required  the  apriori  knowledge  of  the  true  frequency 
values  to  sort  the  signal  frequency  from  the  extraneous  frequencies,  knowledge  that  may 
not  be  available  in  an  automated  method.  The  IV  technique  is  an  iterative  technique  that 
uses  the  past  iteration's  estimation  to  generate  an  improved  set  of  instrumental  variables. 
Instead  of  iterating  to  improve  the  instrumental  variables,  fixed  sets  of  instrumental 
variables  were  used  in  this  example.  The  two  sets  of  instrumental  variables  are  labelled  as 
IVp^rfect  ^^poor  following  tables.  The  IVpgrfg(,{  used  the  exact  signal  as  the 
instrumental  variables  and  represents  the  best  set  of  instrumental  variables  possible.  The 
IVpoor  -'sed  data  generated  with  nine-tenths  of  the  exact  frequency  and  represents  an 

acceptable  set  of  instrumental  variables.  The  other  iterative  algorithms,  MLE  and  GLS, 
were  limited  to  two  and  a  single  iteration,  respectively.  The  number  of  iterations  were 
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limited  because  the  autoregressive  model  of  the  residuals  often  were  identified  as  unstable 
and  the  algorithms  diverged.  The  MPEM  was  overspecified  so  that  there  would  be  four 
extraneous  poles  for  every  signal  pole. 

The  SDOF  data  was  first  corrupted  with  gaussian  random  noise  at  a  signal-to-noise 
ratio  (SNR)  of  40  dB.  Table  2  shows  the  results  from  the  identification  algorithms.  The 
iterative  techniques  (GLS  and  MLE)  failed  to  estimate  the  frequency  and  damping  factors 
locating  the  estimated  poles  on  the  real  axis  in  the  z-plane.  The  instrumental  variables  and 
Pisarenko  methods  estimated  the  damping  factor  as  negative  (unstable).  The  better 
estimates  are  the  overspecified  LS  estimate  and  the  Correlation  Fit.  The  MPEM  method 
estimated  the  damping  factor  as  slightly  negative  which  would  not  have  been  identified  as  a 
signal  pole. 

Next,  correlated  noise  was  added  to  the  signal  in  the  form  of  a  60  Hz  noise  and 
constant  bias.  The  correlated  noise  required  additional  modes  to  model  the  noise.  In  the 
case  of  the  60  Hz  signal  one  addition  mode  is  required.  In  the  case  of  a  constant  term,  half 
a  mode  or  an  increase  of  the  order  by  one  was  needed.  Each  root  of  the  AR  model  is 
equivalent  to  an  exponential  in  the  signal.  An  added  sine  wave  requires  the  addition  of  two 
complex  exponentials;  an  added  constant  requires  a  single  real  exponential.  Table  2  shows 
that  the  MPEM  and  the  overspecified  LS  estimates  are  exact  when  the  model  was 
overspecified. 

An  additional  frequency  was  added  to  the  signal  and  the  same  types  of  noise  were 
also  added.  Tables  3  and  4  show  the  frequency  and  damping  estimates  from  the  seven 
algorithms.  Again  the  best  estimates  were  from  the  overspecified  LS  algorithm  and  the 
MPEM.  Since  the  LS  algorithm  requires  apriori  knowledge  to  distinguish  noise  and  true 
modes  and  the  MPEM  does  it  automatically,  the  MPEM  is  considered  the  better  method. 

The  third  set  of  numerical  experiments  compared  the  MPEM  to  the  LS  and  CF 
methods  when  the  third  frequency  of  25  Hz  was  added.  The  LS  and  CF  methods  represent 
the  best  non-iterative  methods  without  overspecification.  The  overspecified  LS  method 
was  not  considered  because  of  the  inability  to  distinguish  modes  without  apriori 
knowledge.  The  added  noise  was  uncorrelated  gaussian  noise  and  a  higher  signal 
frequency  of  50  Hz  with  an  amplitude  two  orders  of  magnitude  less  than  signal  amplitudes. 
The  higher  frequency  represents  a  higher  mode  whose  amplitude  has  been  attenuated  by  a 
low  pass  filter.  The  final  simulation  contained  the  higher  mode  noise,  uncorrelated 
gaussian  noise  at  40dB,  and  a  constant  bias.  The  final  simulation  represents  real  life  data,  a 
higher  mode,  uncorrelated  noise  at  a  reasonable  level,  and  an  slight  error  in  balancing 
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Table  2:  One  DOF  Simulated  Data 


a. 

II 

o 

X 

N 

^=0.01,  Amp=l.,  (|) 

=10°,  256  pts,  f5= 

2000  Hz) 

Algorithm 

grv  noise 

60Hz 

noise 

constant  bias 

constant  bias 

40dB  rms 

Amp: 

=0.01 

Amp=0.01 

> 

3 

■a 

II 

=0.5 

fd 

fd 

c 

fd 

c 

fd 

c 

LS 

n=2 

* 

0.999 

.012 

0.793 

.054 

n=10 

1.002 

.023 

1.000 

.010 

1.000 

.010 

1.000 

.010 

IV 

poor 

1.012 

-.011 

0.942 

.002 

1.000 

.011 

0.974 

.060 

perfect 

1.008 

-.006 

0.911 

-.023 

0.999 

.012 

0.958 

.084 

PHD 

0.972 

-.002 

1.002 

-.029 

0.999 

.012 

0.793 

.051 

MLE  (2  iter) 

♦ 

* 

♦ 

0.999 

.012 

0.797 

.083 

GLS(  liter) 

* 

* 

* 

★ 

0.999 

.012 

0.793 

.054 

CF 

0.999 

.011 

1.005 

.010 

0.999 

.007 

1.077 

.238 

MPEM  (L=10) 

p=2 

0.994 

-.000* ** 

0.997 

.001 

0.999 

.011 

0.791 

.034 

p=3 

•  t 

.t 

•  t 

.t 

.t 

.t 

1.000 

.010 

_ gr4 

•  t 

•t 

1.000 

.010 

•  t 

•  t 

.t 

JU 

.  1 

*  unable  to  identify 

**  identified  as  extraneous 
t  case  not  considered 
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Table  3:  2D0F  Simulated  Data  Results  for  1st  Mode 
(fd=1.0  Hz,  C=0.01,  Amp=l.,  4>  =10°,  256  pts,  fs=2000  Hz) 


Algorithm  grv  noise  60Hz  noise  constant  bias  constant  bias 


40dB  rms _ Amp=0.01 _ Amp=0.01 _ Amp=0.5 


fd 

c 

fd 

C 

fd 

c 

fd  c 

LS 


n=4 

* 

* 

4c 

.* 

0.999 

.012 

0.791 

.052 

n=10 

1.003 

.014 

1.000 

.010 

1.000 

.010 

1.000 

.010 

IV 

poor 

1.065 

-.612 

0.432 

.915 

1.000 

.011 

0.976 

.057 

perfect 

1.040 

-.595 

1.033 

.672 

0.999 

.012 

0.959 

.083 

PHD 

* 

.* 

4.627 

0.066 

0.999 

.012 

0.790 

.049 

MLE  (2  iter). 

* 

.* 

41 

* 

0.999 

.012 

0.797 

.111 

GLS  (liter) 

* 

* 

4> 

* 

0.999 

.012 

0.791 

.052 

CF 

_  Kc 

* 

* 

0.999 

.007 

1.096 

.217 

MPEM  (L=20) 

p=4 

0.997 

.006 

0.998 

.006 

0.999 

.011 

0.788 

.027 

p=5 

•  t 

.t 

.t 

.t 

•  t 

.t 

1.000 

.010 

_ 

.t 

.t 

1.000 

.010 

•  t 

•  t 

.t 

.t 

*  unable  to  identify 
t  case  not  considered 
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Table  4:  2D0F  Simulated  Data  Results  for  2nd  Mode 
(fd=10.0  Hz,  C=0.01.  Amp=l.,  <])  =20°,  256  pts,  fs=2000  Hz) 


Algorithm 

grv  noise 

40dB  rms 

60Hz  noise 

Amp=0.01 

constant  bias 

Amp=0.01 

constant  bias 

Amp=0.5 

fd 

c 

fd 

C 

fd 

fd 

LS 

n=4 

10.001 

.070 

* 

* 

10.000 

.010 

9.999 

.010 

n=10 

10.000 

.010 

10.000 

.010 

10.000 

.010 

10.000 

.010 

IV 

poor 

10.370 

.020 

9.636 

.051 

10.000 

.010 

9.994 

.010 

perfect 

9.939 

.013 

9.958 

.016 

10.000 

.010 

9.999 

.010 

PHD 

9.904 

.012 

14.969 

0.002 

10.000 

.010 

Q.999 

.010 

MLE  (2  iter) 

8.870 

.108 

9.880 

.172 

10.000 

.010 

10.000 

.010 

GLS  ( 1  iter) 

10.002 

.070 

9.885 

.173 

10.000 

.010 

9.999 

.0 !  0 

CF 

10.002 

.015 

9.997 1 

.0166 

10.000 

.010 

10.000 

.010 

MPEM  (L=20) 

p=4 

9.998 

.010 

10,000 

.010 

10.000 

.010 

10.000 

.010 

p=5 

X 

4. 

.  » 

X 
.  f 

.t 

.t 

X 
.  1 

10.000 

.010 

_ P=^ 

.t 

X 
,  { 

10.000 

.010 

X 
.  1 

.t 

X 

X 

*  unable  to  identify 
t  case  not  considered 
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strain-gage  accelerometers.  Table  5  shows  the  results  from  the  simulated  experiments. 
The  MPEM  estimates  were  quite  accurate. 

2.  Experimental  Free  Response 

The  time  domain  experimental  free  response  from  a  SDOF  system  and  a  cantilever 
beam  cited  by  reference  23  was  used  to  further  compare  the  algorithms.  The  SDOF  system 
was  similar  to  the  base  motion  system  shown  in  Figure  1,  except  the  base  mass  was 
constrained.  The  other  mass  was  constrained  to  a  motion  along  a  single  axis.  The 
damping  was  primarily  provided  by  the  bearing  friction  and  was  definitely  non-viscous. 
The  cantilever  beam  was  made  from  steel  and  was  37.25  inches  long  with  a  cross  section 
of  1.5  by  0.0625  inches.  The  damping  was  structural.  All  of  the  algorithms  except  the 
IV  estimation  were  used. 

The  results  are  shown  in  Table  6.  All  of  the  methods  were  able  to  estimate  the 
natural  frequency  of  the  SDOF  system  reasonably  well,  but  there  were  considerable  errors 
in  the  damping  estimation  when  overspecification  was  not  used.  Log  decrement 
procedures  in  reference  23  estimated  the  damping  factor  to  be  around  2%.  The  MPEM,  the 
overspecified  LS,  the  PHD,  and  CF  damping  estimates  were  all  quite  good. 

Only  the  CF,  the  classical  (non-overspecified)  LS,  and  the  MPEM  were  used  to 
identify  the  cantilever  beams  natural  frequencies  and  damping  factors.  The  LS  method 
failed  to  identify  the  first  r  xle.  The  CF  and  the  MPEM  estimated  the  natural  frequencies 
well,  but  the  MPEM  estimate  of  the  damping  in  the  first  mode  was  lower  than  the  CF 
method  and  on  the  same  order  of  magnitude  as  the  other  modes.  The  MPEM  estimate  is 
probably  the  better  of  the  two  (the  true  values  are  unknown,  since  the  data  is  experimental 
and  a  log  decrement  could  not  be  performed). 

3.  Overspecification  Example 

The  results  from  the  simulated  and  experimental  data  show  the  need  for 
overspecification  and  the  MPEM  method  to  get  an  accurate  estimate  of  the  natural 
frequencies  and  damping  factors  and  hence  an  accurate  estimate  of  the  AR  parameter 
vector.  An  additional  set  of  simulations  were  conducted  using  one  simulated  free  response 
record  to  show  the  effect  of  overspecification.  The  overspecification  size  was  varied  with 
the  LS  algorithm;  and  the  PEM  and  MPEM  algorithms  were  used  to  demonstrate  the  effects 
of  overspecification  and  the  minimum  norm  solution  in  the  z-plane.  Hie  signal  contained 
frequencies  at  10,  20,  and  .30  Hz  all  damped  at  1%.  The  sample  rate  was  chosen  to  be  100 
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Table  5:  3DOF  Simulated  Data 

(fd=l,10,25;  C-01..01..01;  Amp=],l,l;  256  pts,  fs=2000  Hz) 

Algorithm  grv  noise  higher  mode  grv,  4th  mode,  bias 

40  dB  rms  Amp=.01  40  dB,  A=.01,A=.01 

fd 

% 

fd 

LS  (n=6) 

mode  1 

* 

* 

* 

* 

* 

* 

mode  2 

5.242 

.362 

9.974 

.144 

5.207 

.369 

mode  3 

24.741 

.032 

25.066 

.012 

24.729 

.032 

CF 

mode  1 

* 

* 

* 

* 

* 

* 

mode  2 

9.814 

.019 

10.032 

.016 

9.822 

.020 

mode  3 

24.988 

.010 

25.013 

.010 

24.992 

.010 

MPEM  (L=30,  p=6) 

mode  1 

1,001 

.010 

1.000 

.010 

1.000 

.011 

mode  2 

10.000 

.010 

10.000 

.010 

10.000 

.010 

mode  3 

25.001 

.010 

25.000 

.010 

25.001 

.010 

*  unable  to  identify 
t  case  not  considered 
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Table  6:  Experimental  Data  (512  pts,  fs=l(XX)  Hz) 

Algorithm 

SDOF 

1st  mode 

MDOF;  cantilever  beam 

2nd  mode 

3rd  mode 

fd  c 

fd  c 

fd  c 

fd  C 

Analytical 

Estimated 

-1.8  -.020 

1.45 

9.10 

25.49 

LS 

n=2 

1.797 

.156 

.t 

•  t 

.t 

•  t 

•  t 

.t 

n=6 

•  t 

•  t 

* 

* 

8.789 

.003 

24.64 

.003 

n=10 

1.806 

.022 

•  t 

•  t 

.t 

.t 

.t 

.t 

PHD 

1.809 

.023 

.t 

•  t 

•  t 

.t 

.t 

.  T 

MLE  (2  iter) 

1.922 

.008 

.t 

.t 

.t 

•  t 

.t 

.  *r 

GLS 

(liter) 

1.797 

.156 

.t 

•  t 

.t 

.t 

.  T 

CF 

1.804 

.021 

1.515 

.075 

8.805 

.002 

24.66 

.002 

MPEM 

L=10,  p-2 

1.805 

.020 

•  t 

.t 

.t 

.t 

.t 

.t 

L=30,  p=6 

t 

t 

1.512 

.004 

8.804 

.002 

24.66 

.002 

*  unable  to  identify 
t  case  not  considered 
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Hz  (512  points)  so  that  the  system  poles  would  be  well  separated  in  the  z-plane.  Random 
gaussian  noise  was  added  to  make  the  SNR  equal  to  20  dB. 

Figure  26  shows  the  z-plane  with  the  unit  circle  drawn.  The  true  locations  of  the 
poles  are  shown  as  the  circles  in  the  figure  and  the  estimated  pole  locations  from  a  sixth 
order  LS  estimation  are  drawn  as  cross  hairs.  The  estimated  poles  are  very  inaccurate. 
The  estimated  frequencies  are  therefore  inaccurate  as  summarize  in  Table  7.  Figure  27 
shows  what  happens  when  the  order  of  the  LS  estimation  is  increased  to  twelve.  There  are 
estimated  locations  very  close  to  the  true  locations.  Note  that  if  the  true  locations  were  not 
known,  there  would  be  no  way  to  distinguish  the  noise  poles  from  the  signal  poles.  The 
frequency  estimates  are  reasonably  accurate  (1-2.5%  error)  but  the  damping  estimates  are 
inaccurate  (200-400%  error).  The  damping  error  envelopes  given  in  Figures  21,  23,  and 
25  confirm  that  slight  perturbation  in  the  estimated  locations  will  cause  the  large  error  in 
damping  estimation  when  the  damping  is  low  (1%  in  this  case).  The  order  of  the 
estimation  is  increased  further  to  30  to  increase  the  accuracy.  Figure  28  shows  the 
resulting  estimated  pole  locations  and  Table  7  shows  that  both  the  frequency  and  damping 
estimates  are  very  accurate.  Notice  in  Figure  28  that  the  extraneous  roots  are  scattered 
about  the  unit  circle.  When  the  PEM  is  used  with  the  same  order,  the  frequency  and 
damping  estimates  are  about  the  same,  but  Figure  29  shows  that  the  extraneous  roots  aic 
scattered  in  a  more  orderly  fashion.  Still  there  is  no  way  to  distinguish  poles  with  apriori 
knowledge;  even  with  the  "orderly"  position  of  the  poles,  it  would  be  impossible  to 
separate  the  signal  poles.  Figure  30  shows  what  happens  when  the  MPEM  is  used.  The 
true  locations  of  the  poles  are  placed  outside  the  unit  circle,  and  the  estimated  signal  pole 
locations  are  placed  outside  also,  but  the  minimum  norm  solution  places  the  extraneous 
poles  all  inside  the  unit  circle.  Thus,  the  signal  poles  can  be  differentiated  from  the 
extraneous  poles  and  the  high  accuracy  of  oversjjecification  is  achieved. 

4.  Example:  Closely  Spaced  Frequencies 

The  MPEM  allows  for  overspecification,  system  order  determination,  and  automatic 
sorting  of  signal  and  extraneous  poles.  In  practice,  one  would  input  the  free  response  data 
to  a  MPEM  program,  grossly  overspecify  the  system  order,  evaluate  the  singular  values  of 
the  Qm^Qm  niatrix  and  specify  how  many  of  the  principal  eigenvectors  to  use.  The 

MPEM  should  locate  the  signal  poles  outside  the  unit  circle.  The  number  of  principal 
eigenvectors  (PE)  should  be  the  same  number  as  the  true  order  of  the  system,  in  the 
frequency  range  of  interest,  and  the  number  of  signal  poles.  But  if  the  number  of  principal 
eigenvectors  is  slightly  overspecified,  the  MPEM  should  determine  the  true  order  by 
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Figure  26:  Estimated  pole  location  from  the  LS  algorithm  (n=6). 
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Table  7:  Simulated  Data  (512  pts,  fs=200  Hz) 

(fjj  =  10,  20,  30  Hz,  1%  damping  20  dB  grv  noise  added) 


Algorithm _ 1st  mode _ 2nd  mode _ 3rd  mode 

fd  C _ fd  C  fd  c 


LS 


n=6 

12.38 

.227 

28.38 

.092 

83.21 

.128 

n=12 

10.08 

.038 

20.48 

.039 

30.36 

.021 

n=30 

10.00 

.010 

20.01 

.010 

30.00 

.010 

PEM  (L=30,  p=6) 

10.00 

.010 

20.01 

.010 

29.99 

.010 

MPEM  (L=30,p=6) 

10.00 

.010 

20.00 

.010 

29.99 

.010 

75 


Rqq] 


Figure  27;  Estimated  pole  location  from  the  overspecified 
LS  algorithm(n=12). 
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Figure  28:  Estimated  pole  location  from  the  overspecified 
LS  algorithm(n=30). 
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Figure  29:  Estimated  pole  location  from  the  PEM  algorithm  (n=30). 
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Figure  30:  Estimated  pole  location  from  the  MPEM  algorithm  (n=30). 
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locating  the  proper  number  of  poles  outside  the  unit  circle.  If  there  is  too  much  noise,  and 
the  order  has  not  been  overspecified  to  a  large  enough  value,  the  MPEM  method  may  fail. 
Recall  the  example  of  the  two  closely  spaced  frequencies  and  the  singular  value  spectrum 
shown  in  Figure  14.  One  may  chose  either  4,  6,  or  8  as  the  system  order  after  viewing  the 
spectrum.  Also,  the  noise  singular  values  are  very  close  in  magnitude  to  the  signal  singular 
values  and  one  would  expect  that  the  effect  of  the  noise  would  be  great  even  if  the  proper 
order  was  chosen.  Figure  31  shows  the  singular  value  spectra  of  the  Qm^  Qm  matrix  for 

overspecification  orders  of  30  and  50.  The  spectra  shown  in  the  figure  are  intermediate 
results  from  the  MPEM.  Table  8  shows  the  frequency  and  damping  estimates  from  the 
MPEM  when  the  number  of  principal  eigenvectors  are  chosen  as  4,  6,  or  8  for 
overspecification  orders  of  30  and  50.  When  the  overspecifieation  order  is  too  low 
(L=30),  the  MPEM  averages  the  two  closely  spaced  modes  regardless  of  the  selection  of 
the  number  of  principal  eigenvectors  .  The  numerical  difference  between  the  signal 
singular  values  and  the  noise  singular  values  is  low  with  the  low  overspecification  order. 
When  the  order  is  higher  (L=50),  the  selection  of  the  number  of  principal  eigenvectors  does 
have  an  effect  upon  the  estimates.  Notice  that  the  numerical  difference  between  the  signal 
singular  values  and  noise  singular  values  is  higher.  The  closely  spaced  modes  are  averaged 
as  one  value  when  four  principal  eigenvectors  are  used.  The  correct  number  of  modes  are 
identified  when  6  or  8  is  chosen  as  the  number  of  principal  eigenvectors.  Theoreticallv, 
the  number  of  principal  eigenvectors  should  not  be  estimated  too  high,  since  the  mechanism 
for  distinguishing  the  signal  and  extraneous  poles  may  be  lost. 

Overview 

The  key  difficulty  in  estimating  an  AR  model  and  equivalently  estimating  the 
frequencies  and  damping  factors  of  a  system  is  the  effect  of  noise  bias.  In  the  examples, 
the  effect  of  overspecifieation  upon  discrete  time  series  models,  in  this  case  the  AR  model, 
is  shown.  The  difficulty  with  the  overspecified  algorithms  is  the  sorting  of  the  noise 
modes  from  the  true  modes  of  the  model.  The  MPEM  is  able  to  provide  the  benefits  of 
overspecifieation  and  distinguish  between  the  signal  poles  and  the  extraneous  poles.  An 
accurate  low-order  AR  model  can  be  developed  from  the  MPEM  estimates  of  signal  poles. 
However,  the  goal  of  this  report  is  to  identify  the  means  of  constructing  the  whole  discrete 
time  transfer  function.  One  method  to  deal  with  estimation  bias  is  to  simply  identify  the 
overspecified  ARMA  model  in  a  single  stage.  The  alternate  is  a  two-stage  procedure.  The 
first  stage  would  use  free  or  the  impulse  response  and  the  MPEM  to  identify  order  and 
produce  an  accurate  AR  model  containing  signal  frequency  components  only.  The  second 
stage  would  then  have  to  identify  the  MA  parameters  given  the  AR  model  and  the  input- 
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Number 


Figure  31 :  Singular  Value  Spectra  of  the  backward  autocorrelation  matrix 

for  two  closely  spaced  modes. 


Table  8:  MPEM  estimates  for  MDOF  data  with  closely  spaced  modes 


Trial 

L 

P 

mode  1 

mode  2 

mode  3 

freq 

damp 

freq 

damp 

freq 

damp 

# 

Hz 

% 

Hz 

% 

Hz 

% 

1 

4 

14.79 

1.25 

★ 

* 

30.01 

0.88 

2 

30 

6 

14.79 

1.25 

* 

★ 

30.01 

0.88 

3 

8 

14.80 

1.54 

♦ 

* 

29.88 

0.57 

4 

50 

4 

14.79 

1.59 

* 

* 

29.61 

0.54 

5 

50 

6 

14.63 

0.84 

14.98 

0.72 

30.01 

0.90 

6 

8 

14.63 

0.84 

14.98 

0.72 

30.01 

0.91 

True  Values 

14.60 

1.00 

15.00 

1.00 

30.00 

l.Ov) 

*  unable  to  identify 
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output  time  history  set.  The  next  section  explores  the  two  methods  using  simulated  and 
experimental  data. 
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APPLICATIONS 


In  this  section  two  algorithms  used  to  parameterize  a  full  ARMA  model  are 
compared.  The  first  algorithm  is  a  single-stage,  possibly  overspecified,  Least  Squares 
algorithm  and  will  be  denoted  as  the  SLS  (for  single-stage  least-squares)  algorithm.  The 
other  is  a  two  stage  procedure  that  uses  free  response  and  the  MPEM  to  identify  system 
order  and  produce  an  accurate  AR  model  containing  signal  frequency  components  only, 
and  an  MA  LS  procedure  to  identify  the  MA  parameters.  The  two-stage  algorithm  will  be 
denoted  by  2LS.  The  2LS  algorithm  will  produce  a  lower  order  model  without  noise 
modes.  The  SLS  algorithm  will  provide  higher  order  models,  due  to  overspecifration  and 
requires  less  information  since  a  free  response  record  is  not  needed.  The  2LS  algorithm 
requires  additional  data  and  provides  more  information  about  the  system.  It  identifies  the 
signal  frequencies  and  dampings  along  with  providing  a  reduced  order  model  which  will 
allow  for  more  rapid  calculation  of  the  system  response.  Simulated  MDOF  data  and 
experimental  SDOF  and  MDOF  data  will  be  used  to  compare  the  algorithms.  Since  the 
second  stage  of  the  2LS  algorithm  has  not  yet  been  discussed,  it  will  be  presented  here. 

MA  Parameter  Identification 

The  ARMA  model  consists  of  two  parts,  the  Autoregressive  part  and  the  Moving 
Average  part.  While  the  Autoregressive  part  of  the  ARMA  has  an  easily  interpreted 
physical  meaning,  the  Moving  Average  (MA)  part  of  the  model  does  not.  The  MA  pan  of 
the  ARMA  model  relates  the  input  directly  to  the  output.  It  is  the  nonhomogeneous  pan  of 
the  ARMA  model.  A  Moving  Average  process  may  be  modeled  as 

y(k)  =  (bg  -t-  b,z’4  ...+  bpz  ’’)  u(k)  (78) 

A  Moving  Average  process  is  a  weighted  average  of  another  process. 

If  the  AR  parameters  for  a  given  ARMA  model  are  known,  then  the  MA  parameters 
can  then  be  identified  in  a  separate  set  of  calculations  or  "second  stage."  Any  of  ARMA 
identification  algorithms  can  be  modified  to  identify  the  MA  parameters.  But,  for  simplicity 
in  this  repon,  the  MA  parameters  are  determined  by  a  modified  fomi  of  the  LS  algorithm. 
Given  the  input-output  time  histories  and  the  AR  parameters,  the  ARMA  model  can  be 
reduced  to  a  MA  form.  In  order  to  identify  the  MA  parameters,  an  overdetermined  set  of 
simultaneous  equations  is  written  in  matrix  form  as 
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The  LS  solution  of  Eqn  79  is 


=  (82) 

Here  the  notation  of  persistent  excitation  is  clearly  demonstrated.  If  the  scaled  discrete  time 
input  autocorrelation  matrix,  is  singular,  then  the  MA  parameters  can  not  be 

determined.  A  persistent  excitation  provides  full  rank  of  the  input  autocorrelation  matrix 
(reference  16). 

The  second  stage  of  the  2LS  algorithm  solves  the  set  of  equations  through  the  use 
of  a  LS  algorithm.  It  should  be  noted  that  any  of  the  major  algorithms  could  be  used  to 
estimate  the  MA  parameters.  The  MA  LS  parameter  estimates  are  biased,  just  as  the  full 
model  LS  estimate  is  biased.  Other  algorithms  were  considered  for  estimation  of  the  MA 
parameters,  but  since  the  LS  algorithm  is  the  simplest  of  the  major  algorithms  to 
implement,  it  was  chosen  for  this  sample  study. 

Examples 

A  series  of  experiments  were  petfonned  to  illustrate  the  application  of  the  proposed 
two-stage  approach.  Simulated  data  for  a  MDOF  system  was  used  to  compare  the 
algorithms.  Also,  data  for  a  SDOF  and  a  MDOF  system  were  acquired  experimentally 
from  actual  structural  systems  and  used  to  examine  the  2LS  algorithm.  Results  show  that  a 
model  developed  using  single-stage  LS  parameter  estimates  (without  overspecification) 
fails  to  predict  response  adequately  for  both  the  SDOF  and  MDOF  system,  while  the  model 
developed  using  ihe  2LS  algorithm  parameter  estimates  provides  very  good  response 
predictions.  Overspecified  models  can  be  identified  using  the  SLS  algorithm  that  can  be 
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used  to  accurately  predict  response.  All  measured  quantities  cited  in  the  following 
discussion  are  accelerations  measured  in  g's. 

1 .  Simulated  Data 

Input  and  output  time  histories  were  simulated  to  compare  the  two  candidate 
algorithms,  the  2LS  and  the  SLS.  The  input  was  a  pseudo-random  binary  input  as  shown 
in  Figure  32.  A  sixtli  order  arbitrary  ARMA  model  was  used  to  generate  the  system  output 
time  history  for  this  input.  The  output  was  corrupted  with  gaussian  random  noise  at  a 
signal-to-noise  ratio  of  40  dB.  The  corrupted  response  is  shown  in  Figure  33.  The  SLS 
algorithm  was  used  to  identify  a  sixth  and  twelfth  order  ARMA  model.  The  FRF 
magnitudes  of  the  identified  and  exact  ARMA  models  are  shown  in  Figure  34.  The  Figure 
illustrates  that  even  a  modest  level  of  noise  effects  the  SLS.  The  the  first  mode  of  the  exact 
transfer  function  is  accurately  modeled  by  the  twelfth  order  model  and  second  mode  is 
modeled  but  not  as  accurately.  The  sixth  order  model,  recall  the  actual  order  of  the  system 
is  six,  is  not  represented  accurately.  Notice  the  strange  behavior  in  the  FRF  magnitude  plot 
of  the  12th  order  model  past  30  Hz.  The  extraneous  modes  are  the  cause  of  this  behavior. 

The  noise  used  in  the  example  has  a  flat  spectral  density  over  the  entire  frequency 
range.  If  these  data  were  experimentally  acquired,  low  pass  filters  would  have  limited  the 
bandwidth  of  the  signal  and  the  bandwidth  of  the  noise.  The  equivalent  digital  filter  of  a 
four  pole  low  pass  Butterworth  filter  (reference  5)  was  used  to  simulate  the  effect  of  actual 
low  pass  filters.  The  frequency  response  of  the  digital  filter  is  shown  in  Figure  35.  If  the 
input  and  output  are  filtered  by  similar  low  pass  filters,  the  transfer  function  (within  the 
bandwidth)  will  be  unaffected  by  the  filters,  but  the  filters  will  truncate  the  high  frequency 
noise.  Figures  36  and  37  are  the  simulated  filtered  input  and  output  produced  when  the 
input  and  output  time  histories  are  modified  by  the  digital  filter  (actually  an  ARMA  model). 
When  the  SLS  algorithm  is  used  with  the  filtered  input-output  data  sets,  the  results  are 
much  better.  Figure  38  shows  the  FRF  magnitude  plots  of  the  6th  and  12th  order  models. 
Both  models  follow  the  first  two  modes  closely,  but  the  third  mode  is  not  identified. 

The  FRF  of  the  2LS  models  identified  using  the  filtered  and  unfiltered  data  are 
shown  in  Figure  39.  The  2LS  models  are  6th  order  and  were  obtained  by  assuming  that 
the  AR  part  of  the  model  is  known  precisely.  The  model  obtained  from  the  unfiltered  data 
is  inaccurate.  The  model  ubiained  trom  the  tiuerea  Oaia  is  accurate  for  the  first  two  modes 
and  third  mode  is  identified  but  not  very  accurately.  Certain  characteristics  of  the  two 
algorithms  are  readily  illustrated  by  the  example.  First,  the  SLS  models  improve  as  the 
order  is  increa.sed,  and  secondly,  the  algorithms  work  best  when  the  noise  is  band-limited. 
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Figure  32: 


Pseudo-random  binary  input. 
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Figure  33:  Noisy  response  to  the  pseudo  random  binary 

input. 
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Figure  34;  Frequency  Response  Function  magnitude 
for  SLS  models  (unfiltered  data  used). 
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Figure  35;  Frequency  Response  of  digital  filter. 
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Figure  36;  Filtered  pseudo-random  binary  input. 
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Figure  38:  Frequency  Response  Function  magnitude 
for  SLS  models  (filtered  data  used). 


If  the  AR  portion  of  the  model  can  be  identified  accurately,  the  2LS  algorithm  can  produce 
the  lowest  possible  order  model  which  is  as  accurate  as  a  much  higher  order  SLS  model. 
The  overspecified  SLS  models  will  have  dynamics  past  the  bandwidth  of  interest.  A  high 
order  SLS  and  a  2LS  model  can  be  compared  on  a  FRF  magnitude  plot  to  compare 
accuracy. 

2.  The  SDOF  experiment 

A  mechanical  system  similar  to  that  depicted  in  Figure  1  was  used  to  acquire  input- 
output  time  histories  to  test  the  2LS  identification  procedure.  The  system  was  fabricated 
using  two  masses  which  were  constrained  to  one-dimensional  movement  by  linear  bearings 
and  two  parallel,  cylindrical  rods.  The  masses  were  connected  by  springs.  The  left-hand 
or  base  mass  could  be  moved  manually  to  provide  input  to  the  system;  the  response  of  the 
free  mass  was  considered  the  output  of  the  system.  The  experiment  corresponds  to  a 
SDOF  base  motion  problem.  The  system  damping  was  provided  by  the  bearing  friction 
and  is  nonviscous. 

The  2LS  procedure  requires  the  free  response  for  AR  parameter  estimation  in  the 
first  stage.  Experimental  data  was  acquired  by  constraining  the  base  mass  and  displacing 
the  free  mass  from  equilibrium.  A  strain-gage  accelerometer  was  used  to  measure  the  low 
frequency  acceleration  response.  A  typical  free  response  is  shown  in  Figure  40.  This 
record  contains  500  data  points  with  a  sample  frequency  of  100  Hz,  so  there  is  a  total  of  5 
seconds  of  free  response  data.  The  ordinary  LS  estimates  for  the  damped  natural  frequency 
and  damping  factor  are  1.820  Hz  and  23%  for  this  data  set.  The  MPEM  estim.ates  (L=10) 
are  1 .806  Hz  and  2.0%.  The  frequency  estimates  are  similar  but  the  damping  factors  differ 
by  an  order  of  magnitude.  A  log  decrement  procedure  was  used  to  estimate  the  damping 
factor  (between  1.5  and  1.8%),  verifying  that  the  MPEM  estimate  of  the  damping  factor  is 
more  reasonable.  Other  experimental  free  response  records  were  acquired  and  repeatability 
of  the  MPEM  estimates  verified.  The  average  damped  natural  frequency  and  damping 
factor  from  four  MPEM  estimates  were  1.805  Hz  and  2.1%.  The  average  estimates  were 
used  to  produce  the  AR  parameters  used  in  the  second  stage.  Experimental  input-output 
time  histories  were  obtained  by  manual  excitation  of  the  base  mass.  The  input  to  the 
system  was  measured  as  the  acceleration  of  the  base  mass  by  a  second  strain-gage 
accelerometer.  Figures  41  to  43  show  three  different  inputs  (500  pts,  fs=100  Hz),  Case  A, 

Case  B,  and  Case  C,  respectively,  which  were  applied  to  the  system.  The  single-stage  LS 
(SLS)  model  estimates  (without  overspecification)  obtained  from  the  three  input-output 
time  histories  were  inconsistent,  and  the  resulting  ARMA  models  provided  inaccurate 


timo  (see) 


Figure  40:  Typical  free  response  of  the  SDOF  system. 
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Figure  42:  Base  acceleration,  Case  B. 
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Figure  43;  Base  acceleration.  Case  C. 


99 


predictions  of  system  response.  The  predicted  responses  are  obtained  through  the  use  of 
the  estimated  ARMA  models  and  the  input  time  histories  without  knowledge  of  the  actual 
output  time  histories.  Figure  44  shows  the  predicted  response  (in  comparison  to  the  actual 
response)  to  the  input  Case  B  using  the  SLS  model  obtained  from  input-output  time 
histories  for  Case  A.  The  predicted  response  is  typical  of  the  cases  considered  with 
models  developed  using  the  SLS  estimates.  Figures  45  to  47  show  the  predicted  responses 
(in  comparison  to  the  actual  responses)  to  input  Cases  A,  B,  and  C,  respectively,  from  the 
2LS  model  obtained  from  the  input-output  time  histories  for  Case  A.  The  2LS  model 
obtained  using  the  input-output  time  history  corresponding  to  input  Cases  B  and  C  provide 
similar  predicted  responses.  The  magnitude  frequency  response  plots  for  the  three  2LS 
ARMA  models  are  shown  in  Figure  48  while  the  frequency  response  plots  for  the  three 
SLS  ARMA  models  are  shown  in  Figure  49.  The  SLS  ARMA  models  have  grossly 
overestimated  the  damping.  The  2LS  ARMA  models  are  very  similar  to  each  other  when 
compared  using  the  FRF,  and  appear  to  be  more  reasonable  than  the  SLS  models. 

Overspecification  of  the  model  in  the  SLS  procedure  improves  the  damping 
estimates  but  introduces  extraneous  poles  and  zeroes  as  shown  in  Figure  50.  The  damping 
and  frequency  estimates  are  given  in  Table  9  as  the  model  order  is  increased  in  the  SLS 
procedure.  The  damped  natural  frequency  is  about  4  percent  of  the  nyquist  frequency,  and 
one  would  expect  high  errors  in  the  damping  factor  estimates.  The  response  prediction  of 
input  Case  B  using  the  tenth  order  SLS  ARMA  model  obtained  from  input  Case  A  is 
shown  in  Figure  51.  Note  that  the  damping  estimate  is  still  too  high  and  there  is  some  error 
in  frequency. 

3.  The  MDOF  Experiment 

A  cantilevered  thin  steel  beam  (length=37.25",  width=1.5",  thickness=  0.0625") 
was  used  as  a  MDOF  experimental  system.  The  beam  was  cantilevered  on  a  block,  the 
block  was  fixed,  and  the  beam  was  initially  deformed  to  provide  free  vibration.  A  strain- 
gage  accelerometer  was  mounted  on  the  tip  of  the  beam  to  measure  free  response.  The 
analog  data  was  filtered  (with  cutoff  frequencies  of  16,  31.5,  and  63  Hz)  using  an  analog 
four-pole  low-pass  Butterworth  filter  and  then  digitized.  The  MPEM  algorithm  was  used 
to  identify  the  natural  frequencies  and  damping  factors  from  the  free  vibration  response. 
The  results  appear  in  Table  10  along  with  the  analytical  frequency  values  computed  for  a 
cantilevered  beam.  Note  that  the  founh  mode  was  identified  in  Trial  5  as  an  extraneous 
mode.  In  this  case,  the  noise  perturbed  the  pole  locations  (corresponding  to  the  founh 
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Figure  44:  Second  order  SLS  model  prediction 
of  SDOF  response  to  input  Case  B. 
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Figure  45:  Second  order  2LS  model  prediction 
of  SDOF  response  to  input  Case  A. 
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Figure  46:  Second  order  2LS  model  prediction 
of  SDOF  response  to  input  Case  B. 


103 


Rasponsc  (g) 


Figure  47:  Second  order  2LS  model  prediction 
of  SDOF  response  to  input  Case  C. 
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Figure  48:  Magnitude  frequency  response  plots 
of  second  order  2LS  models. 
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Figure  49:  Magnitude  frequency  response  plots 
of  second  order  SLS  models. 
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Figure  50:  Magnitude  frequency  response  plots 
of  overspecified  SLS  models  (input  Case  A). 
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Table  9:  Frequency  and  damping  estimates 
from  single  stage  LS  estimates  of  Case  A 

Order 

Freq 

c 

(Hz) 

(%) 

2 

1.600 

13.2 

4 

1.800 

5.6 

6 

1.798 

3.9 

8 

1.798 

3.8 

10 

1.793 

3.1 

12 

1.794 

3.1 

20 

1.788 

2.8 
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Figure  51:  Tenth  order  2LS  model  prediction 
of  SDOF  response  to  input  Case  B. 
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Table  10:  MPEM  estimates  for  MDOF  exp.  data 


Trial 

L 

mode  1 

mode  2 

mode  3 

mode  4 

freq 

damp 

ffeq 

damp 

freq 

damp 

freq  damp 

# 

Hz 

% 

Hz 

% 

Hz 

% 

Hz  % 

1 

30 

1.516 

.500 

8.801 

.223 

(filtered) 

(filtered) 

2 

30 

1.518 

.594 

8.801 

.229 

(filtered) 

(filtered) 

3 

40 

1.512 

.322 

8.804 

.205 

24.66 

.168 

(filtered) 

4 

40 

1.512 

.454 

8.800 

.271 

24.64 

.181 

(filtered) 

5 

50 

1.517 

.686 

8.793 

.307 

24.57 

.201 

(unidentified) 

6 

50 

1.518 

.463 

8.800 

.274 

24.55 

.330 

48.29  .034 

AVG 

1.515 

.503 

8.800  .251 

24.60 

.220 

48.29 

.034 

Analytical 

1.45 

9.10 

25.49 

49.95 

mode)  inside  the  unit  circle,  indicating  that  more  trials  are  necessary.  The  average  of  the 
identified  values  were  then  used  in  the  second  stage  of  the  2LS  procedure. 

The  base  block  that  the  beam  was  cantilevered  from  was  attached  to  a  shaker  to 
provide  controllable  system  input.  The  input  was  measured  as  the  acceleration  response  of 
the  base  block  using  a  strain  gage  accelerometer  and  the  system  response  was  measured  as 
the  beam  tip  acceleration.  The  voltage  sent  to  the  shaker  was  a  random  binary  voltage 
generated  by  computer.  The  voltage  was  randomly  either  +1V  or  -  IV  for  a  discrete  time. 
Two  sets  of  filtered  input  accelerations  to  this  system  are  shown  in  Figures  52  and  53 
(Case  D  and  Case  E,  500pts,  fs=250  Hz),  Most  of  the  input  energy  was  contained 
between  the  frequencies  of  approximately  5  and  35  Hz.  The  shaker  was  unable  to  provide 
input  at  frequencies  as  low  as  the  first  mode,  and  therefore  the  first  mode  was  not  excited. 
The  MA  parameters  for  this  system  were  estimated  from  the  data  assuming  either  two  or 
three  significant  modes. 

The  frequency  response  plots  of  SLS  ARMA  models  of  order  4,  6,  8,  and  10 
obtained  using  data  from  input  Case  D  are  shown  in  Figure  54.  As  the  overspecification  is 
increased,  the  peaks  are  better  defined.  2LS  models  of  order  four  and  six  were  also 
obtained  from  input  Case  D  using  the  frequency  and  damping  estimates  of  just  the  second 
and  third  mode  for  the  fourth  order  model  and  the  estimates  of  the  second,  third  and  fourth 
modes  for  the  sixth  order  model  to  determine  the  AR  parameters.  The  frequency  response 
plots  for  these  two  cases  are  shown  in  Figure  55.  The  frequency  response  plot  shows  that 
the  sixth  order  model  is  approximately  equivalent  to  the  fourth  order  model  in  the  frequency 
range  of  interest.  The  response  predictions  from  input  Case  E  using  the  sixth  and  tenth 
order  SLS  models  are  shown  in  Figures  56  and  57  in  comparison  with  the  actual  system 
response.  The  sixth  order  model  does  not  predict  response  as  well  as  the  overspecified 
model.  The  fourth  order  2LS  model  predicts  the  response  as  well  as  the  overspecified  SLS 
model  as  shown  in  Figure  58,  and  the  prediction  can  be  developed  in  approximately  40% 
of  the  time  required  for  the  tenth  order  overspecified  model. 

Overview 

The  comparison  of  results  from  the  numerical  simulations  show  the  improved 
transient  response  predictions  from  using  models  developed  with  the  2LS  method.  The 
results  from  experimental  systems  are  quite  good  even  though  the  experimental  system  is 
nonviscously  damped,  and  the  ARMA  model  is  better  suited  for  a  viscously  damped 
system.  The  MPEM  was  successful  in  identifying  the  frequency  and  damping  factors  of 
the  first  four  modes  for  the  cantilever  beam.  The  MPEM  did  fail  to  provide  estimates  of  the 
fourth  mode  of  the  experimental  MDOF  system  in  Trial  5.  In  this  case,  the  measured  free 
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Figure  52:  Acceleration  input  to  the  MDOF  system,  input  Case  D. 


112 


Input  (g) 


Figure  53:  Acceleration  input  to  the  MDOF  system,  input  Case  E. 
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Figure  54;  Magnitude  frequency  response  plots 
of  SLS  models  (input  Case  D). 
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Figure  55;  Magnitude  frequency  response  plots 
of  2LS  models  (input  Case  D). 
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Figure  56.  Sixth  order  SLS  model  used  to  predict  response 
of  the  MDOF  system  to  input  Case  E. 
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Figure  57.  Tenth  order  SLS  model  used  to  predict  response 
of  the  MDOF  system  to  input  Case  E. 
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Figure  58:  Fourth  order  2LS  model  used  to  predict  response 
of  the  MDOF  system  to  Case  E. 
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response  contained  too  much  noise.  The  MPEM  sensitivity  toward  noise  can  be 
diminished  by  further  overspecification  of  the  model  but  that  was  not  attempted  in  this 
study.  An  overspecified  ordinary  LS  estimate  of  the  AR  parameters  could  have  been  used 
to  identify  the  frequencies  and  dampings,  but  it  would  have  been  difficult  to  distinguish  the 
signal  roots  from  the  extraneous  roots. 

The  deficiencies  of  an  ordinary  single-stage  identification  scheme  arise  from  bias 
problems  in  the  estimation  of  the  AR  parameters  of  the  ARMA  model.  The  advantage  of  the 
two-stage  procedure  is  in  the  determination  the  AR  parameters  separately  from  the  MA 
parameters  allowing  greater  accuracy  in  the  estimation.  The  2LS  method  reduces  the  bias 
in  the  AR  parameters  before  the  second  stage  and  allows  for  the  identification  of  a  reduced 
order  ARMA  model.  The  disadvantage  of  the  2LS  method  is  that  it  requires  an  additional 
set  of  data,  a  free  response  record.  An  alternative  is  to  form  the  impulse  response  function 
using  time  domain  data,  but  this  usually  requires  an  ensemble  average  of  multiple  tests. 
Ideally,  the  single-stage  algorithm  can  be  used  with  overspecified  systems  and  the 
extraneous  poles  and  roots  can  be  eliminated  providing  a  small  order  accurate  model.  The 
difficulty  is  in  the  automatic  sorting  of  the  extraneous  poles  and  roots  from  the  signal  pole 
and  roots. 
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NONLINEAR  MODELS 


Discrete  time  series  can  be  used  for  transient  response  prediction  of  linear 
structures.  Free  vibration  can  be  modeled  by  a  special  form  called  an  Autoregressive  (AR) 
model.  When  structural  nonlinearities  are  present,  it  may  be  possible  to  modify  the  form  of 
the  AR  model  to  incorporate  the  nonlinearities.  One  possibility  is  to  allow  the  model 
parameters  to  become  functions  of  state.  This  report  explores  the  proper  form  of  the 
parameter  functions  for  various  nonlinear  structures.  Two  numerical  case  studies  and 
experimental  results  are  used  to  evaluate  the  model  form.  This  section  presents  the  results 
of  a  preliminary  study  intended  to  investigate  the  suitability  of  modeling  nonlinear  structural 
systems  with  Discrete  Time  Series  models. 

The  Nonlinear  AR  Model 

Free  vibration  for  linear  structures  can  be  modeled  as  an  Autoregressive  process. 
The  Autoregressive  (AR)  models  can  be  used  to  extract  modal  parameters  such  as  damping 
factors  and  natural  frequencies.  The  AR  models  can  be  parameterized  from  measured  free 
vibration  using  the  identification  algorithm  presented  in  the  previous  sections.  The  models 
can  then  be  used,  if  initial  conditions  are  known,  to  predict  vibration.  The  standard  AR 
model  for  a  single  output  system  is 

x(k)  =  -aj  x(k-l)  -  a2x(k-2)  -...-apX(k-p)  +  e(k)  (83) 

where  x(k)  is  the  system's  free  response  at  discrete  time  k,  the  a's  are  the  parameters  of  the 
model,  e(k)  is  the  residual  or  error,  and  p  is  the  system  order.  For  multiple  output 
systems,  x  becomes  a  vector  of  responses  and  the  parameters  become  matrices.  The  AR 
model  can  be  put  into  the  form 

[1  +  ajZ’'  +  a2Z'^+...+apZ’’’]  x(k)  =e(k)  (84) 

using  the  definition  of  the  forward  shift  operator,  z.  The  AR  model  used  to  predict  the  free 
response  of  a  linear  system  can  be  modified  in  an  attempt  to  incorporate  certain  types  of 
system  nonlinearities.  In  this  report  the  constant  coefficients  of  the  AR  model  become 
functions  of  state  in  order  to  form  the  nonlinear  Autoregressive  (NLAR)  models.  The 
general  NLAR  model  is 

[1+  a  j(x)  z  +  ..+  a  p(x)  z  x(k)  =e(k)  (85) 
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where  the  x  vector  represents  the  state 


X  =  [x(k-l),x(k-2),...,x(k-p)]  (80 

The  NLAR  model  is  an  approximation.  The  accuracy  of  the  NLAR  model  depends  upon 
the  magnitude  of  nonlinearity,  the  form  of  individual  AR  functionals,  the  accuracy  in  the 
parameterization  of  the  functionals,  and  the  initial  conditions.  The  form  of  the  NLAR 
model  has  to  be  determined  before  the  predictive  model  is  parameterized.  Single-degree-of- 
freedom  (SDOF)  nonlinear  oscillators  were  studied  to  gain  insight  into  model  forms.  The 
general  NLAR  model  for  a  nonlinear  SDOF  oscillator  is 

[1+  aj(x)  z  V  a  2(x)  z  x(k)  =  e(k)  (87) 

The  studies  focused  upon  the  form  of  the  functionals  aj  and  a2. 

Examples 

Two  case  different  systems  were  studied  using  simulated  time  series  data.  These 
systems  were  a  Duffing  Oscillator  and  a  viscous-coulomb  damped  oscillator.  Noiseless 
position  time  history  data  were  generated  numerically  by  integrating  the  equations  of 
motions  via  the  fourth  order  Runge-Kutta  algorithm. 

System  response  data  for  a  complex  nonlinear  mechanical  system  was  also  acquired 
experimentally.  A  Piper  Apache  oleo-pneumatic  landing  gear  strut  was  used.  Experimental 
position  data  were  acquired  using  the  oleo-pneumatic  strut  test  facility  shown  in  Figure  59. 
Regardless  of  its  source,  the  position  data  were  used  to  estimate  the  parameters  of  the 
proposed  models  in  a  least  squares  (LS)  fashion.  The  first  two  points  from  the  position 
data  were  then  used  along  with  the  models  to  generated  predicted  response.  The  position 
data  and  the  predicted  response  were  compared  to  qualitatively  compare  the  models.  The 
investigations  were  focused  upon  the  evaluation  of  the  model  form  and  not  necessarily 
upon  the  identification  algorithms  or  the  influence  of  noise.  For  these  reasons  the 
simulation  data  were  not  corrupted  with  noise,  but  noise  in  the  experimental  data  was 
unavoidable. 

l.The  Duffing  Oscillator 

The  Duffing  Oscillator  represents  an  SDOF  system  where  the  stiffness  is  nonlinear 
and  the  damping  is  vi.scous.  The  equation  of  motion  for  a  Duffing  Oscillator  is 
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Figure  59:  Oleo-pneumatic  shock  strut  experimental  apparatus. 


x  +  cx  +  ax  +  px^  =  0 


(88) 


If  p  is  positive,  the  oscillator  becomes  more  stiff  as  the  displacement,  x,  increases  and  the 
oscillator  is  said  to  have  a  hard  spring.  If  P  is  negative,  the  oscillator  becomes  less  stiff  as 
the  displacement  increases,  and  the  oscillator  is  said  to  have  a  soft  spring.  Only  the  results 
using  hard  springs  will  be  presented  in  this  report.  NLAR  models  for  the  Duffing  oscillator 
have  been  cited  in  the  literamre  in  the  past.  Reference  24  suggests  that  the  NLAR  model 

x(k)  +  (bj+Cj  exp( -x^(k-l)))  x(k-l)  +  aj  x(k-2)  =  E(k)  (89) 

could  be  used  to  approximate  the  Duffing  Oscillator.  This  model  will  be  referred  to  as  the 
exponential  model,  since  reference  24  has  incorporated  the  stiffness  nonlinearity  of  the 
Duffing  oscillator  into  an  exponential  functional  for  aj.  Usually  differential  equations  can 

be  approximated  through  finite  differencing  procedures.  Here,  central  differencing  of  Eqn 
88  is  used  to  suggest  an  alternate  and  somewhat  simpler  nonlinear  form, 

x(k)+  (  bj  +  C|  x2(k-l))  x(k-l)  +  a2(k-2)  =  e(k)  (90) 

This  NLAR  model  will  be  referred  to  as  the  parabolic  model,  since  ai  has  a  parabolic 
functional  form.  Also,  the  linear  AR  model,  where  the  parameter  are  constants,  can  also  be 
used  to  approximate  the  Duffing  Oscillator  and  the  influence  of  the  nonlinearity 
investigated. 

The  numerical  simulations  used  the  following  form  of  the  hard  spring  Duffing 
Oscillator 

X  +  0.25  X  +  14.8  X  +34.02  x^  =0  (91) 

which  is  also  the  oscillator  used  by  reference  24  .  The  initial  conditions  and  time  increment 
used  in  the  numerical  integration  were 

x(0)=l.  X  (0)  =0.  At  =0.01  sec  (92) 

The  three  models,  the  exponential  NLAR  model,  the  parabolic  NLAR  model,  and  the  AR 
model,  were  estimated  from  500  points  of  data.  The  linear  model,  identified  by  the 
Modified  Principal  Eigenvectors  Method,  the  exponential  model,  identified  by  a  LS 
algorithm,  and  the  parabolic  model,  identified  by  a  LS  algorithm,  are  all  given  in  Table  1 1 . 
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Table  11:  Models  used  for  Duffing  Oscillator 


Trial  AR  Model  Exponential  Model  Parabolic  Model 

_ ai _ 32 _ ai _ bi _ a2 _ ai  bj  a2 


1  -1.9946  0.9978  -1.9910  -0.00545  0.9975  -1.9960  0.00339  0.9975 

2  _ I _ * _ -1.9862  -0.0153  0.9975  -1.9960  0.00338  0.9975 

*  not  considered 
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Figure  60  shows  the  predictions  of  the  three  models  (using  the  first  two  simulated  data 
points  as  initial  conditions)  in  comparison  to  the  simulated  data.  Both  NLAR  models 
pyerform  bettci  than  the  lineai'  model.  The  differences  between  the  models  and  the  analytical 
Duffing  Oscillator  can  be  examined  by  plotting  the  damped  natural  frequency  versus 
displacement  x,  as  shown  in  Figure  61.  The  percent  of  critical  damping  as  a  function  of 
displacement  can  also  be  examined,  but  is  not  included  in  this  report.  The  natural 
frequency  of  the  Duffing  Oscillator  is  (a  +  Px2)l/2.  The  natural  frequency  of  the  NLAR 
models  are  determined  by  linearizing  the  models  with  x(k-l)  equal  to  x.  Figure  61  shows 
that  the  exponential  model  is  a  good  approximation  for  .2<  I  x  I  <  1 ,  while  the  parabolic 
model  is  good  for  much  higher  amplitudes. 

Simulated  data  were  used  to  test  if  the  performance  of  the  NLAR  models  degrade  as 
the  amplitude  is  increased  (essentially,  increasing  the  nonlinearity).  The  equation  of  motion 
was  integrated  again  with  the  initial  position  value  doubled.  The  models  were  then  used 
with  the  new  initial  conditions  to  predict  the  response.  Figure  62  shows  the  response 
predictions  of  the  NLAR  models  in  comparison  to  simulated  data.  The  exponential  model 
fails  to  predict  the  nonlinear  response.  The  parabolic  model's  response  prediction  is  good 
within  the  resolution  of  the  plot.  The  results  suggest  that  the  exponential  model  does  not 
approximate  the  amplitude  dependence  of  the  Duffing  Oscillator  as  well  as  the  parabolic 
model.  This  can  be  confirmed  by  estimating  the  parameters  of  the  models  using  the  new 
simulated  Duffing  Oscillator  data.  The  form  of  an  NLAR  model  should  remove  the 
amplitude  dependence  of  the  parameters.  The  NLAR  models  estimated  from  the  new  data 
are  also  given  in  Table  11.  The  exponential  model's  parameters  change  while  the  parabolic 
model  i.s^variant  to  the  change  in  initial  conditions. 

2.  Viscous-Coulomb  Damped  Oscillator 

The  numerical  studies  also  included  a  viscous-coulomb  damped  oscillator.  The 
equation  of  motion  for  the  viscous-coulomb  damped  oscillator  is 

X -I- 2^0)  X  +  co^x c^  sign(x)  =0  (93) 

where  c^  is  the  coulomb  damping  coefficient.  No  NLAR  model  for  this  system  was  found 
in  the  literature.  Finite  differencing  of  Eqn  93  suggests  an  NLAR  model  of  the  form 
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Figure  60:  Duffing  oscillator  response  and  predictions 
from  NLAR  and  AR  models. 
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Figure  61 :  Frequency  of  the  Duffing  oscillator  and  the  models 
as  a  function  of  displacement. 
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Figure  62:  Duffing  oscillator  response  and  response  predictions 
from  NLAR  models  (larger  initial  amplitude). 
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x(k)  +  Ej  x(k-l)  +  a2x(k-2)  +bjSign  [x(k-l)  -  x(k-2)]  =  e(k) 


(94) 


The  linear  AR  mcxlel  may  also  be  used  if  the  coulomb  damping  term  is  insignificant.  The 
viscous-coulomb  damped  oscillator  used  in  the  simulation  was 

X  -1-  0.5  X  +  25  X  +  0.5  sign(x)  =0  (95) 

The  initial  conditions  and  time  increment  used  in  the  numerical  integration  were  the  same  as 
Eqn  92.  The  linear  model,  identified  by  the  MPEM,  and  the  NEAR  model  identified  by  a 
LS  algorithm  are  given  in  Table  12.  Figure  63  shows  the  predictions  of  the  two  models  in 
comparison  to  the  simulated  data.  Both  predictions  are  very  good  with  the  NEAR  model 
prediction  following  the  simulated  data  more  closely.  The  effectiveness  of  the  linear  model 
in  the  vibration  prediction  shows  the  nonlinearity  of  Eqn  93  may  not  be  great  enough  to 
warrant  an  NEAR  model  for  this  particular  set  of  parameters. 

3.A  Eanding  Gear 

Nonlinear  data  from  a  Piper  Apache  landing  gear  strut  (see  reference  25)  was  used 
to  evaluate  the  NEAR  models.  The  experimental  data  (500  pts,  fs=200  Hz)  appear  in 

Figure  64.  Notice  the  constant  value  of  the  data  past  about  1.5  second.  The  landing  gear 
sticks  at  a  new  equilibrium  after  every  oscillation  due  to  the  strong  effect  of  the  bearing 
friction  in  the  gear.  This  bearing  friction  can  be  modeled  as  coulomb  damping.  The  NEAR 
used  to  model  the  vibration  was  of  the  form 

x(k)  +  Ej  x(k-l)  +  a2x(k-2)  +  b,sign  [x(k-l)  -  x(k-2)]+b2  =  e(k)  (96) 

which  is  the  same  as  the  viscous-coulomb  damped  NEAR  model  plus  a  constant  term  to 
represent  the  nonzero  equilibrium  position.  An  AR  model  has  no  mechanism  to  simulate 
coulomb  sticking,  so  only  the  NEAR  model  was  identified  by  a  ES  algorithm.  The  NEAR 
model  identified  at  the  200  Hz  sample  rate  is  not  accurate,  possibly  due  to  noise  bias. 
However,  when  the  data  is  decimated  by  using  only  every  fifth  data  point  (reducing  the 
sample  rate  to  40  Hz),  the  effect  of  the  noise  bias  is  not  as  great.  The  resulting  NEAR 
model  is  given  in  Table  12,  and  the  response  prediction  of  the  NEAR  model  is  shown  in 
Figure  64.  The  NEAR  prediction  is  more  accurate  than  the  predictions  of  the  oscillation 
given  by  reference  25  and  exhibits  the  sticking  phenomenon  associated  with  coulomb 
damping. 
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Table  12:  Mcxiels  used  for  Viscous-Coulomb  Damped  Oscillator 


a2  bi _ b2 


SinmlatLoa 

AR  model 

NLAR  model 

-1.9904 

-1.9924 

0.9929 

0.9949 

4.77E-5 

Experimental 

Strut 

-1.9156 

0.9512 

0.00175 

-0.08677 
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Response  (in) 


Figure  63:  Viscous-coulomb  damped  oscillator  response 
and  predictions  from  NLAR  and  AR  models. 
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Disp] acement  (in) 


tlma  (sec) 


Figure  64:  Experimental  data  and  NEAR 
and  AR  model  response  predictions. 
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Overview 

Finite  differencing  techniques  were  used  to  suggest  forms  for  nonlinear 
autoregressive  models.  The  NLAR  models  when  parameterized  by  LS  techniques  using 
noiseless  simulated  data  provided  excellent  prediction.  Linear  AR  models  provided 
response  predictions  nearly  as  good  as  the  NLAR  models  when  the  degree  of  the 
nonlinearity  was  not  great.  The  exponential  NLAR  model  is  accurate  when  the  degree  of 
nonlinearity  of  the  Duffing  Oscillator  is  small.  But,  the  parabolic  model  is  a  better  model 
for  the  Duffing  Oscillator  particularly  for  large  amplitude  response.  Experimental  data  from 
an  oleo-pneumatic  strut  was  also  used  to  evaluate  NLAR  models.  The  prediction  of  the 
experimental  response  by  the  NLAR  was  very  reasonable  and  is  better  than  the  predictions 
by  the  nonlinear  differential  equations  parameterized  by  reference  25. 

The  use  of  time  domain  discrete  models  in  the  form  of  Autoregressive  Moving 
Average  (ARMA)  models  has  been  shown  to  be  advantageous  in  the  response  prediction  of 
forced  linear  vibration.  When  the  structure  is  nonlinear,  the  ARMA  models  may  have  to  be 
extended  to  model  the  nonlinearities  just  as  the  AR  models  are  extended  to  NLAR  models. 
The  new  nonlinear  ARMA  model's  parameters  will  be  functions  of  state.  This  is  a  topic 
requiring  further  research. 
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CONCLUSIONS  AND  RECOMMENDATIONS 


Discrete  time  series  models  for  both  linear  and  nonlinear  systems  have  been  shown 
to  be  an  effective  and  accurate  means  of  representation  and  prediction.  Linear  models  are 
effective  when  the  order  is  large  enough  and  the  sample  rate  is  at  least  an  order  of 
magnitude  higher  than  the  upper  bandwidth  limit.  The  time  series  models  become  more 
accurate  as  the  sample  rate  increases,  but  the  effect  of  noise  perturbations  on  the  parameters 
also  become  greater  as  the  sample  rate  increases.  Some  trade-off  must  be  made  when 
selecting  a  sample  rate  for  acquiring  data  to  be  used  in  parameter  identification.  The  sample 
rate  must  be  high  enough  to  insure  accurate  predictions  once  the  model  is  parameterized. 
But  the  sample  rate  must  be  low  enough,  so  that  the  model  can  be  parameterized  in  the 
presence  of  noise. 

There  are  various  techniques  to  reduce  the  effect  of  noise  bias  ranging  ftx>m  iterative 
techniques  to  overspecification  techniques.  The  iterative  techniques  often  experience 
convergence  problems  when  an  interim  model  is  identified  as  unstable.  The 
overspecification  techniques  are  easier  to  implement  but  result  in  excessively  large  model 
order.  The  model  must  be  reduced  to  the  proper  order  by  some  technique  if  it  is  to  be 
effective  in  response  prediction  or  control.  Time  domain  modal  methods  compare  mode 
shapes  using  modal  confidence  factors  and  are  able  to  eliminate  the  noise  model  from  the 
system  model  in  the  identified  model.  The  difficulties  with  modal  methods  are  they  require 
multiple  measurement  locations,  are  time  consuming,  and  are  not  well  suited  for  automated 
data  processing.  Ideally,  one  would  like  to  measure  a  limited  number  of  inputs  and  outputs 
to  parameterize  the  model.  When  a  SISO  ARMA  model  is  desired,  only  the  single  input 
and  the  single  output  time  histories  should  be  needed  to  characterized  the  model.  The 
traditional  signal  processing  techniques  utilize  single  input  and  single  output  time  histories 
but  provide  no  means  of  separating  noise  from  signal  modes  in  the  identified  model.  One 
has  to  use  the  entire  model,  with  noise  and  signal  fractions,  to  predict  vibration  as  was 
done  with  the  overspecified  SLS  models.  The  2LS  algorithm  allows  the  AR  part  of  the 
model  to  be  reduced  to  its  proper  size  through  the  use  of  the  MPEM.  A  non -overspecified 
LS  estimation  of  the  MA  parameters  completes  the  algorithm.  The  MPEM  was  shown  to 
be  a  highly  accurate  and  automatic  method  for  determining  system  order.  The  MPEM 
identifies  natural  frequencies  and  eliminates  noise  bias  from  free  response  records  even  in 
the  case  of  closely  spaced  modes.  The  main  drawback  of  the  2LS  algorithm  is  that  it 
requires  the  additional  free  response  record.  The  frequency  and  damping  error  envelopes 
demonstrate  the  high  accuracy  required  when  identifying  the  autoregressive  part  of  the 
model. 


Experimental  data  from  a  SDOF  and  a  MDOF  was  obtained,  and  ARMA  models 
were  successfully  identified  using  the  SLS  and  the  2LS  algorithms.  The  2LS  models,  with 
orders  of  the  proper  size,  successfully  predicted  forced  dynamic  response.  The  SLS 
models  had  to  be  much  larger,  greatly  overspecified,  to  be  as  accurate  as  the  2LS  models. 
The  FRF  plots  graphically  show  how  the  SLS  models  begin  to  follow  the  first  few  modes 
of  vibration  and  how  all  modes  are  accurately  represented  when  the  order  is  overspecified 
to  large  values.  The  FRFs  also  show  the  dynamics  of  the  noise  modes  beyond  the 
bandwidth  limit.  There  are  also  noise  modes  in  the  bandwidth,  but  their  effects  cannot  be 
uncovered  because  their  amplitudes  are  small  when  compared  to  the  signal  modes.  There 
does  not  seem  to  be  any  ill-effects  of  carrying  the  noise  modes  in  the  model  when  it  is  used 
for  vibration  prediction,  other  than  the  additional  computational  time  required.  But,  it  is 
conceivable  that  there  are  some  contributions  in  the  predicted  response  from  these  noise 
modes.  The  2Lo  models  have  no  noise  modes  and  may  be  used  effectively  for  response 
prediction  purposes,  and  also  for  control  algorithm  purposes  where  pole  and  zero 
cancellations  are  required. 

The  optimum  algorithm  would  require  only  the  single  input  single  output  time 
histories  and  would  overspecify  the  model  and  eliminate  the  extraneous  noise  modes.  A 
two- stage  algorithm  would  first  have  to  identify  the  AR  model  from  an  overspecified 
ARMA  model  and  then  identify  the  MA  parameters  in  a  LS  fashion  without  requiring  the 
addition  set  of  free  response  data.  The  other  drawback  of  the  2LS  method,  as  it  is 
proposed,  is  that  the  MA  pan  of  the  model  is  not  overspecified.  There  is  some  noise  in  the 
input  and  overspecification  would  allow  the  MA  parameters  to  be  better  identified.  But, 
there  is  no  way  of  identifying  which  part  of  the  MA  model  is  extraneous  without  the  use  of 
the  AR  parameters.  The  best  method  would  identify  an  overspecified  ARMA  model  and 
split  it  into  a  modal  form  using  partial  fraction  expansion,  and  identify  the  true  modes  of  the 
system  using  some  criteria.  The  problem  is  what  type  of  criteria  to  use. 

Real  nonlinear  free  response  was  used  to  evaluate  the  NLAR  models.  The  NLAR 
model  predicts  the  free  vibration  of  an  oleo-pneumatic  strut  fairly  well.  However,  an 
accurate  model  requires  the  careful  selection  of  the  data  sampling  rate  and  the  determination 
of  the  influence  of  measurement  noise  on  the  parameters.  This  has  not  been  considered  in 
this  investigation.  The  algorithm  that  identified  the  NLAR  models  was  a  non-overspecified 
LS  algorithm  which  produces  biased  parameter  estimates.  The  bias  at  the  lower  sample  rate 
had  a  smaller  effect  than  the  bias  at  the  high  sample  rate.  No  algorithms  were  investigated 
to  reduce  the  bias.  Given  further  research  in  the  area,  the  NLAR  models  could  become 
easier  to  parameterize  and  more  attractive  to  use. 
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APPENDIX  A  -  Alternative  Formulation 

The  governing  differential  equation  of  any  linear  system  can  be  put  into  first  order  state 
space  form 

x(t)  =  A  x(t)  +  B  u(t)  (Al) 

y(t)  =  Cx(t)  (A2) 

The  differential  equation,  Eqn  Al  can  be  rearranged  for  solution  after  multiplying  by  exp(- 
At)  to 

^[e'^^^ft)]  =e’^*B  u(t)  (A3) 

The  solution  of  Eqn  A3  is  made  by  simply  integrating  between  tj^  and  t 

t 

e'^‘  x(t)  -  e  x(tj^)  =J  B  u(x)  dx  (A4) 

•k 

After  rearranging  and  multiplying  by  exp(At),  the  solution  is 

t 

x(t)  =  e'^^^’^^xCtj^)  +  J  B  u(x)  dx  (A5) 

S. 

Since  the  selection  of  t  and  tj(  is  arbitrary,  let  t|j  be  kh  and  t  be  (k+l)h.  Notice  that  the  input 
can  vary  between  the  sample  times  t  and  tj^  and  that  Eqn  A5  is  the  exact  solution  but  the 
integral  depends  upon  the  input  time  function.  However,  if  one  does  not  know  the  input 
over  the  entire  time  domain,  but  only  at  discrete  time  instance,  Eqn  A5  can  be 
approximated.  If  we  assume  the  input  is  constant  over  the  interval  between  t  and  t^  and  is 
equal  to  uCtj^),  then  Eqn  A5  can  be  reduced  to 

h 

x([k-*-l]h)  =  e'^‘'  x(kh)+  [j  e^^  dt  ]  Bu(kh)  (A6) 

o 

or 
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x([k+l  ]  h)  =  <D  x(kh)  +  r  u(kh)  (A7) 

Equation  A7  allows  one  to  find  an  approximate  solution  of  Eqn  A1  at  discrete  time 
instances  by  performing  a  recursive  matrix  operation.  The  matrix  operation  can  be  put  into 
a  more  compact  form  through  the  use  of  the  forward  shift  operator  as 

(z  I  -  O)  x(kh)  =  r  u(kh)  (A8) 

Finally  the  solution  is  found  through  inversion  and  the  use  of  Eqn  A2  as 

y(kh)  =  C  (zI.<D)"  V  u(kh)  (A9) 

here,  the  discrete  time  transfer  function  matrix  is 

G(z)=C(zI.(D)‘^  r  (A  10) 
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APPENDIX  B:  The  AR  model 


Complex  Exponentials 


Consider  a  summation  of  exponentials 

x(t)  =  exp  t) 

k=l 


(Bl) 


where  and  can  be  complex  or  real.  If  x(t)  represents  the  response  of  a  structure, 
should  appear  as  complex  conjugates  pairs  and  the  corresponding  has  to  appear  as 
complex  conjugates.  If  the  sample  interval  is  h,  then  the  response  at  j  sample  instances 
later  is 


x(t+jh)  =  ^  b^  exp  (X^  (t+jh)) 


which  can  be  written  as 


k=l 


x(t+jh)  =  z’x(t)  ^  b^  exp(Xj^t) 
k=l 


(B2) 


(B3) 


where  z  is  forward  shift  operator  and 

zk=exp(Xich) 

Writing  Eqn  B3  for  j=0,l,2,...  m  in  matrix  form 


■*  X(t) 

-1 

-1 

...  -1 

zx(t) 

-^1 

• 

m  .  . 

m 

m 

m 

_z  x(t) 

r  1  A 

bjexp(X|t) 


V 


V’‘P<V> 


=  0 


J 


or 


(B4) 


(B5) 


'P  P  =  0.  (B6) 

For  non-trivial  b^,  the  determinant  of  'F  is  zero.  The  determinant  in  terms  of  the  first 
column  of  'F  is 


zm  x(t)  Cm  +  x(t)  Cm-i  + . +  z  x(t)  Ci  +  x(t)  Cq  =0  (B7) 


147 


where  Cj  is  co-factor  of  the  i+1  element  in  column  1  of 'F.  Dividing  through  by  Cj^  and 
by  one  can  write 


{ 1  +  aj  z'l  +  a2  z‘2  +....  +  )  x(t)  =0  (B8) 

Equation  B8  is  the  equation  for  an  Autoregressive  process,  the  polynomial  in  z  is  the 
characteristic  equation  in  the  time  domain.  It  is  apparent  from  Eqn  B8  that  a  summation  of 
exponentials  (real  or  complex)  is  an  Autoregressive  process.  Also,  for  Eqn  B8  to  hold, 
'F  must  be  singular.  If  Xj  ^  Xj  (that  is  zj  5^  zj  ),  'F  becomes  singular  when  the  first 
column  is  a  multiple  of  any  other  column;  this  occurs  when  z  e  {  Z2^  2^  )  .  The 

roots  of  the  characteristic  equation  in  the  time  domain  are  given  by  Eqn  B4  which  can  be 
mapped  to  the  roots  of  the  characteristic  equation  in  the  Laplace  domain  and  eventually  to 
the  modal  frequencies  and  damping  factors. 

Alternative  Formulation 

The  equations  of  motion  for  any  linear  homogeneous  equation  can  be  put  into  state 
space  form  as 

i(t)  =  A  x(t)  (B9) 

y(i)  =  Cx(t)  (BIO) 

where  x(t)  is  the  vector  of  states.  The  exact  solution  of  Eqn  B8  is  given  by  Eqn  A8  when 
r  is  zero 

(zl  -  O)  x(kh)  =  0  (Bll) 

for  non-trivial  x,  the  determinant  of  (zl  -  0)  must  be  zero.  The  determinant  results  in  the 
characteristic  equation  in  the  Time  Domain.  Thus,  every  state  must  follow  an 
Autoregressive  process  and  since  y(t)  is  a  linear  combination  of  the  states  at  a  given  time, 
y(t)  is  also  an  Autoregressive  process. 
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APPENDIX  C:  Singular  Value  Decomposition 


Any  matrix.  A,  belonging  to  can  be  decomposed  into  a  product  of  three  matrices, 

(Cl) 

(C2) 

and  where  the  Oj  are  termed  the  singular  values  of  A.  The  singular  values  are  ordered  from 
the  largest  to  the  smallest.  The  columns  of  U  are  denoted  by  uj  and  the  columns  of  V  are 
denoted  by  vj.  The  matrices  U  and  V  are  orthogonal  to  themselves,  that  is 

u‘u  =  U  U  =  I 

V  v;=  I 

V  V  =  I  (C3) 

If  A  is  of  rank  r,  then 

o,>Oj>...>o,>o„,=  ...=  o„-0  (C4) 

the  nullspace  of  A  is 

N(A)  =  span  ....vJ  (C5) 

and  the  range  of  A  is 

R(A)  =  span  [Uj,  u^. ....  uj  (C6) 

The  pseudo- inverse  of  A  is 

A#=  (ATa)-UT  =  V  I  -  luT  (C7) 

A  symmetric  matrix  B,  can  be  represented  by 

B  =  a”^A  =  U  I  (V  L  U”^)  =  U  Z  (C8) 

note  that  Z  2  is  also  a  diagonal  matrix.  The  condition  number  of  A  is  the  ratio  of  the 
maximum  singular  value  and  the  minimum  singular  value.  If  the  singular  values  past  the 
rth  singular  value  are  nearly  zero,  then  the  minimum  norm  approximation  of  A  is 


A  =  U  Z  VT 


where 


,nxn 


U  e 
V  €  C 
Z  =diag(a,,a2,...,an) 
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(C9) 


A  =  ^  a.  u.  v.^ 

i=i 

References  12  and  26  provide  complete  discussions  of  Singular  Value  Decomposition 
(SVD). 
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APPENDIX  D-Singular  Value  Decomposition  Solution  of  Normal  Equations 

The  ill-conditioning  problems  associated  with  the  solution  of  an  overdetermined  set 
of  equations  has  been  well  documented.  The  LS  algorithm  requires  the  solution  of  an 
overdetermined  set  of  equations  by  forming  the  Normal  Equations.  Reference  27  states, 
"the  Normal  Equations  are  ill-conditioned  especially  when  the  number  of  parameters  is 
large,"  and  suggests  the  solution  of  the  overdetermined  set  of  equations  by  repeated 
Householder  transformations  until  the  equations  are  in  triangular  form.  Reference  28 
recognizes  the  ill-conditioning  problem  and  proposes  solution  through  singular  value 
decomposition  (SVD)  or  QR  decomposition.  Reference  29  states,  "utilization  of  the  normal 
equations  was  numerically  unstable,  due  to  the  squaring  of  the  condition  number  of  the 
coefficient  matrix,"  and  recommends  the  use  of  the  QR  decomposition  in  the  solution  of  the 
overdetermined  set  of  equations.  The  solution  of  the  overdetermined  set  of  equations  via 
SVD  is  shown  to  be  equivalent  to  solution  of  the  normal  equations  by  reference  12. 

Solution  by  SVD  has  been  shown  to  have  many  useful  properties  including  '•ank 
determination  and  determination  of  the  conditioning  number  of  a  matrix.  However,  the 
SVD  solution  of  the  overdetermined  set  of  equations  will  be  shown  to  have  another  useful 
property,  a  well  conditioned  solution.  Given  the  overdetermined  set  of  equations 


Ae=b 


(HI) 


where  A  is  Nxp.  The  LS  solution  involves  the  normal  equations 

0  =  (ATA)-lAfb 


(D2) 


A  can  be  decomposed  to 


A  =  U  Z  VT  (D3) 

where  U  is  an  Nxp  unitary  matrix,  S  is  a  diagonal  matrix,  and  V  is  a  pxp  orthogonal 
matrix.  Substitution  of  the  decomposition  of  A  into  Eqn  D2  yields 

e  =  (V  I  U  Z  V^)'‘V  Z  b  (D4) 


since  U  is  unitary 
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0  =  (  V  L  ^V^)' V  Z  b  (D5) 

Note  that  VZ  is  the  SVD  of  A'Ta  and  has  the  condition  number  of 

2 

y  ^Inax 

k(A  A)= — Y-  (D6) 

<^mln 

in  terms  of  the  singular  values  of  A.  The  conditioning  of  the  normal  equations  is  the 
square  of  the  conditioning  of  the  overdetermined  equations.  The  inverse  in  Eqn  D5  is  easily 
computed  to  be  VZ'^V^  reducing  Eqn  D5  to 

e=VZ*u'^b  (D7) 

Equation  D7  is  the  solution  to  the  overdetermined  set  of  equation  via  the  pseudo-inverse 
calculated  through  the  use  of  the  SVD  of  A.  Equation  D7  is  the  solution  to  the  Least 
Squares  problem  without  the  formation  of  A^A.  The  condition  number  of  the  pseudo¬ 
inverse  is  the  square  root  of  Eqn  D6.  If  1/jc(AT'A)  approaches  the  precision  of  the 
computer,  the  solution  via  Eqn  D2  will  be  ill-conditioned  while  solution  via  Eqn  D7  may 
not.  If  I/kCA’Ta)  does  not  approach  the  precision  of  the  computer,  Eqn  D2  and  Eqn  D7 
will  produce  the  same  results.  For  example,  if  a  minicomputer  has  15  and  1/2  digits  of 
resolution  in  double  precision  and  k(aT’A)  approaches  10^6  the  solution  of  Eqn  D2  will  be 
ill-conditioned.  The  LS  solution  of  normal  equations  can  sometimes  be  very  ill- 
conditioned.  Solution  of  the  overdetermined  equation  through  the  use  of  the  SVD  can 
alleviate  the  problem.  But  solution  via  SVD  is  costly  in  CPU  time  and  memory  storage. 
In  many  cases  when  the  normal  equation  are  solved  by  Eqn  D2,  A^A  (mxm)  and  A^b 
(mxl)  can  be  formed  directly  form  data  without  first  forming  A  (nxm)  or  b  (nxl)  where  n 
is  the  number  of  overdetermined  equations  and  m  is  the  number  of  parameters.  If  the 
normal  equations  are  solved  by  SVD  the  matrix  A  has  to  be  formed.  Solution  by  (D2) 
requires  inversion  of  a  mxm  matrix.  Solution  by  SVD  requires  the  decomposition  of  a  nxm 
matrix.  If  n  is  many  times  greater  than  m,  then  the  SVD  solution  will  be  considerably 
slower  than  the  solution  via  inversion  of  the  normal  equations. 
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APPENDIX  E:  The  frequency  error  envelopes 


Once  the  roots  of  the  time  domain  characteristic  equation  are  found,  the  frequency 
and  damping  factors  can  be  extracted  from  Eqn  16 


s  =  Co)  +  j  0^  =  |;|-[ln  Izl  +  arg(z)] 


(El) 


The  error  in  the  estimated  damped  natural  frequency  is 

5(qi='J- arg(^j -z)  (E2) 

where  is  the  estimated  pole  location  and  z  is  the  real  location.  The  maximum  error  is 

e  1  e  1  e 

otm  =T-arctanrT-=T-arctan -  (E3) 

^  h  IzT  K  exp(-C£0h) 

where  e  is  the  radius  of  uncertainty.  The  percentage  error  is  found  by  dividing  by  toj  and 
can  be  put  into  non-dimensional  form  by  noting 


(Oj 


(E4) 


The  damped  natural  frequency  percentage  error  function  is 


error  =  ■ 


■arctan 


which  for  small  e  is  linear  in  e. 


-TtC  ^ 

E/exp( — 7^  - ) 


V 


0)^ 


(E5) 
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